diff --git a/.gitignore b/.gitignore
index b669c7f..6c62dc2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -161,6 +161,7 @@ examples/output/
 # IDE
 .vscode/
 .idea/
+.claude/
 *.swp
 *.swo
 *~
@@ -172,4 +173,5 @@ examples/output/
 .Spotlight-V100
 .Trashes
 ehthumbs.db
-Thumbs.db
\ No newline at end of file
+Thumbs.db
+.early.coverage
\ No newline at end of file
diff --git a/MIGRATION_PLAN.md b/MIGRATION_PLAN.md
deleted file mode 100644
index 0d9fcf6..0000000
--- a/MIGRATION_PLAN.md
+++ /dev/null
@@ -1,758 +0,0 @@
-# CFD-Python Migration Plan
-
-This document outlines the required changes to update cfd-python bindings to work with CFD library v0.1.6.
-
-## Current State
-
-- **cfd-python version:** 0.2.0 (compatible with CFD v0.1.6)
-- **Target CFD library:** v0.1.6
-- **Status:** ✅ COMPLETE - All migration phases finished
-
-## What's New in v0.1.6
-
-CFD library v0.1.6 introduces **Modular Backend Libraries** - a major architectural change that splits the library into separate per-backend components:
-
-- `cfd_core` - Grid, memory, I/O, utilities (base library)
-- `cfd_scalar` - Scalar CPU solvers (baseline implementation)
-- `cfd_simd` - AVX2/NEON optimized solvers
-- `cfd_omp` - OpenMP parallelized solvers (with stubs when OpenMP unavailable)
-- `cfd_cuda` - CUDA GPU solvers (conditional compilation)
-- `cfd_api` - Dispatcher layer and high-level API (links all backends)
-- `cfd_library` - Unified library (all backends, backward compatible)
-
-**Key Changes:**
-- Backend availability API for runtime detection
-- `ns_solver_backend_t` enum with `SCALAR`, `SIMD`, `OMP`, `CUDA` values
-- New functions: `cfd_backend_is_available()`, `cfd_backend_get_name()`, `cfd_registry_list_by_backend()`, `cfd_solver_create_checked()`
-- Improved error codes: `CFD_ERROR_UNSUPPORTED` for unavailable backends
-
-**Impact on cfd-python:**
-- No breaking API changes (maintains v0.1.5 compatibility)
-- Can now query available backends at runtime
-- Better error messages when backends unavailable
-- Same linking approach (use `cfd_library` or `CFD::Library` CMake target)
-
-## Breaking Changes Summary
-
-### 1. Type Name Changes
-
-All type names have been changed to follow C naming conventions:
-
-| Old (cfd-python) | New (CFD v0.1.6) |
-|------------------|------------------|
-| `FlowField` | `flow_field` |
-| `Grid` | `grid` |
-| `SolverParams` | `ns_solver_params_t` |
-| `SolverStats` | `ns_solver_stats_t` |
-| `Solver` | `ns_solver_t` |
-| `SolverStatus` | `cfd_status_t` |
-| `SolverCapabilities` | `ns_solver_capabilities_t` |
-| `OutputRegistry` | `output_registry` |
-| `SimulationData` | `simulation_data` |
-
-### 2. Solver Registry API Changes
-
-**Old (Global Singleton):**
-```c
-void solver_registry_init(void);
-Solver* solver_create(const char* type_name);
-int solver_registry_list(const char** names, int max_count);
-```
-
-**New (Context-Bound):**
-```c
-ns_solver_registry_t* cfd_registry_create(void);
-void cfd_registry_destroy(ns_solver_registry_t* registry);
-void cfd_registry_register_defaults(ns_solver_registry_t* registry);
-ns_solver_t* cfd_solver_create(ns_solver_registry_t* registry, const char* type_name);
-int cfd_registry_list(ns_solver_registry_t* registry, const char** names, int max_count);
-```
-
-### 3. Error Handling Changes
-
-**Old:**
-```c
-typedef enum { SOLVER_STATUS_OK, SOLVER_STATUS_ERROR, ... } SolverStatus;
-```
-
-**New:**
-```c
-typedef enum { CFD_SUCCESS, CFD_ERROR_NOMEM, CFD_ERROR_INVALID, ... } cfd_status_t;
-
-// New error functions:
-void cfd_set_error(cfd_status_t status, const char* message);
-const char* cfd_get_last_error(void);
-cfd_status_t cfd_get_last_status(void);
-const char* cfd_get_error_string(cfd_status_t status);
-void cfd_clear_error(void);
-```
-
-### 4. Simulation Data Structure Changes
-
-**New fields added:**
-```c
-typedef struct {
-    // ... existing fields with renamed types ...
-    ns_solver_registry_t* registry;  // NEW: context-bound registry
-    char output_base_dir[512];       // NEW: output directory
-} simulation_data;
-```
-
-### 5. Output Enum Changes
-
-```c
-// Old
-OUTPUT_PRESSURE
-
-// New
-OUTPUT_VELOCITY_MAGNITUDE
-```
-
-### 6. New Solver Types
-
-New solvers added in v0.1.5+ (inherited by v0.1.6):
-
-- `"explicit_euler_omp"` - OpenMP parallel explicit Euler
-- `"projection_omp"` - OpenMP parallel projection
-- `"conjugate_gradient"` - CG linear solver (internal)
-
-### 7. Backend Availability API (v0.1.6)
-
-New in v0.1.6 for runtime backend detection:
-
-```c
-// Backend enum
-typedef enum {
-    NS_SOLVER_BACKEND_SCALAR = 0,
-    NS_SOLVER_BACKEND_SIMD = 1,
-    NS_SOLVER_BACKEND_OMP = 2,
-    NS_SOLVER_BACKEND_CUDA = 3
-} ns_solver_backend_t;
-
-// Backend availability functions
-bool cfd_backend_is_available(ns_solver_backend_t backend);
-const char* cfd_backend_get_name(ns_solver_backend_t backend);
-int cfd_registry_list_by_backend(ns_solver_registry_t* registry,
-                                  ns_solver_backend_t backend,
-                                  const char** names,
-                                  int max_count);
-ns_solver_t* cfd_solver_create_checked(ns_solver_registry_t* registry,
-                                        const char* type_name);
-```
-
----
-
-## Missing Features
-
-### Boundary Conditions (Completely Missing)
-
-The CFD library has a comprehensive BC API that is not exposed:
-
-**BC Types:**
-- `BC_TYPE_PERIODIC` - Periodic boundaries
-- `BC_TYPE_NEUMANN` - Zero-gradient
-- `BC_TYPE_DIRICHLET` - Fixed value
-- `BC_TYPE_NOSLIP` - No-slip walls
-- `BC_TYPE_INLET` - Inlet velocity
-- `BC_TYPE_OUTLET` - Outlet conditions
-
-**Functions to expose:**
-```c
-// Core BC functions
-cfd_status_t bc_apply_periodic(flow_field* field, bc_edge_t edge);
-cfd_status_t bc_apply_neumann(flow_field* field, bc_edge_t edge);
-cfd_status_t bc_apply_dirichlet_scalar(double* field, ...);
-cfd_status_t bc_apply_dirichlet_velocity(flow_field* field, ...);
-cfd_status_t bc_apply_noslip(flow_field* field, bc_edge_t edge);
-cfd_status_t bc_apply_inlet(flow_field* field, bc_edge_t edge, const bc_inlet_config_t* config);
-cfd_status_t bc_apply_outlet_scalar(double* field, ...);
-cfd_status_t bc_apply_outlet_velocity(flow_field* field, ...);
-
-// Backend control
-cfd_status_t bc_set_backend(bc_backend_t backend);
-bc_backend_t bc_get_backend(void);
-bool bc_backend_available(bc_backend_t backend);
-
-// Inlet configuration helpers
-bc_inlet_config_t bc_inlet_config_uniform(double u, double v);
-bc_inlet_config_t bc_inlet_config_parabolic(double u_max, double v_max);
-bc_inlet_config_t bc_inlet_config_custom(bc_velocity_profile_func profile, void* user_data);
-
-// Outlet configuration helpers
-bc_outlet_config_t bc_outlet_config_zero_gradient(void);
-bc_outlet_config_t bc_outlet_config_convective(double advection_velocity);
-```
-
-### Derived Fields & Statistics (Missing)
-
-```c
-// Field statistics structure
-typedef struct {
-    double min_val;
-    double max_val;
-    double avg_val;
-    double sum_val;
-} field_stats;
-
-// Derived fields container
-derived_fields* derived_fields_create(size_t nx, size_t ny);
-void derived_fields_destroy(derived_fields* derived);
-void derived_fields_compute_velocity_magnitude(derived_fields* derived, const flow_field* field);
-void derived_fields_compute_statistics(derived_fields* derived, const flow_field* field);
-```
-
-### CPU Features Detection (Missing)
-
-```c
-typedef enum {
-    CPU_FEATURE_NONE = 0,
-    CPU_FEATURE_SSE2 = (1 << 0),
-    CPU_FEATURE_AVX = (1 << 1),
-    CPU_FEATURE_AVX2 = (1 << 2),
-    CPU_FEATURE_NEON = (1 << 3),
-} cpu_features_t;
-
-cpu_features_t cfd_get_cpu_features(void);
-```
-
----
-
-## Migration Plan
-
-### Phase 1: Fix Breaking Changes (Critical) ✅ COMPLETED
-
-**Priority:** P0 - Must complete before any other work
-
-**Status:** Completed on 2025-12-26
-
-**Tasks:**
-
-- [x] **1.1 Update bundled headers or remove them**
-  - Chose Option A: Removed bundled headers, require installed CFD library
-  - Deleted `src/cfd_lib/` directory
-
-- [x] **1.2 Update type definitions in cfd_python.c**
-  - Replaced all `FlowField` → `flow_field`
-  - Replaced all `Grid` → `grid`
-  - Replaced all `SolverParams` → `ns_solver_params_t`
-  - Replaced all `SolverStats` → `ns_solver_stats_t`
-  - Replaced all `Solver` → `ns_solver_t`
-  - Replaced all `SimulationData` → `simulation_data`
-
-- [x] **1.3 Update solver registry code**
-  - Added module-level `g_registry` handle
-  - Updated to use `cfd_registry_create()` and `cfd_registry_register_defaults()`
-  - Updated `cfd_solver_create()` and `cfd_registry_list()` calls
-
-- [x] **1.4 Update error handling**
-  - Replaced with `cfd_status_t` error codes
-  - Added error handling API: `get_last_error()`, `get_last_status()`, `get_error_string()`, `clear_error()`
-  - Exposed CFD_SUCCESS, CFD_ERROR, CFD_ERROR_* constants to Python
-
-- [x] **1.5 Update simulation API calls**
-  - Updated function signatures to match new API
-  - Uses `derived_fields` for velocity magnitude computation
-
-- [x] **1.6 Fix output enum**
-  - Replaced `OUTPUT_PRESSURE` with `OUTPUT_VELOCITY_MAGNITUDE`
-
-- [x] **1.7 Update CMakeLists.txt**
-  - Added CFD library version check (require >= 0.1.6)
-  - Added `CFD_BUILD_INCLUDE_DIR` for generated export header
-  - Find CFD library headers in correct paths
-  - Link modular backend libraries (cfd_api, cfd_core, cfd_scalar, cfd_simd, cfd_omp, cfd_cuda)
-  - Added GNU linker groups on Linux for circular dependency resolution
-  - Automatic CUDA library detection for optional GPU support
-
-**Actual effort:** 1 day
-
-### Phase 2: Add Boundary Condition Bindings (Important) ✅ COMPLETED
-
-**Priority:** P1 - Required for useful Python API
-
-**Status:** Completed on 2025-12-26
-
-**Tasks:**
-
-- [x] **2.1 Create BC type enums for Python**
-  - Added BC_TYPE_* constants (PERIODIC, NEUMANN, DIRICHLET, NOSLIP, INLET, OUTLET)
-  - Added BC_EDGE_* constants (LEFT, RIGHT, BOTTOM, TOP)
-  - Added BC_BACKEND_* constants (AUTO, SCALAR, OMP, SIMD, CUDA)
-
-- [x] **2.2 Implement core BC wrapper functions**
-  - `bc_apply_scalar(field, nx, ny, bc_type)` - Apply BC to scalar field
-  - `bc_apply_velocity(u, v, nx, ny, bc_type)` - Apply BC to velocity fields
-  - `bc_apply_dirichlet(field, nx, ny, left, right, bottom, top)` - Fixed value BC
-  - `bc_apply_noslip(u, v, nx, ny)` - Zero velocity at walls
-
-- [x] **2.3 Implement inlet BC wrappers**
-  - `bc_apply_inlet_uniform(u, v, nx, ny, u_inlet, v_inlet, edge)` - Uniform inlet
-  - `bc_apply_inlet_parabolic(u, v, nx, ny, max_velocity, edge)` - Parabolic profile
-
-- [x] **2.4 Implement outlet BC wrappers**
-  - `bc_apply_outlet_scalar(field, nx, ny, edge)` - Zero-gradient outlet
-  - `bc_apply_outlet_velocity(u, v, nx, ny, edge)` - Zero-gradient outlet
-
-- [x] **2.5 Implement backend control**
-  - `bc_set_backend(backend)` - Set active backend
-  - `bc_get_backend()` - Get current backend
-  - `bc_get_backend_name()` - Get backend name string
-  - `bc_backend_available(backend)` - Check availability
-
-- [x] **2.6 Add BC tests**
-  - Tested all BC types (Neumann, Dirichlet, no-slip)
-  - Verified backend detection (OMP available, SIMD detected)
-  - All tests pass
-
-**Actual effort:** 1 day
-
-### Phase 2.5: CI/Build System for v0.1.6 (Critical) ✅ COMPLETED
-
-**Priority:** P0 - Required for v0.1.6 compatibility
-
-**Status:** Completed on 2025-12-29
-
-**Tasks:**
-
-- [x] **2.5.1 Implement dual-variant wheel builds**
-  - Matrix build strategy for CPU-only and CUDA-enabled wheels
-  - CPU wheels: Linux, macOS, Windows (Scalar + SIMD + OpenMP backends)
-  - CUDA wheels: Linux, Windows (All CPU backends + CUDA, Turing+ GPUs)
-  - Artifact naming: `wheel-{os}-{variant}` for differentiation
-
-- [x] **2.5.2 Fix CMakeLists.txt for modular libraries**
-  - Link all modular CFD libraries individually
-  - GNU linker groups on Linux for circular dependencies
-  - Automatic CUDA library detection
-
-- [x] **2.5.3 Update CI test infrastructure**
-  - Install CUDA runtime (12.4.0) for CUDA wheel tests
-  - Use standard `pip` instead of `uv` for stable ABI wheel installation
-  - Test matrix: Python 3.9 and 3.13 on all platforms
-  - Use apt-based CUDA installation on Linux (more reliable than runfile)
-
-- [x] **2.5.4 Ensure PEP 427 compliance**
-  - Standard wheel filenames (no variant suffixes)
-  - Variant differentiation through artifact names only
-  - PyPI-compatible wheel naming
-
-**Actual effort:** 1 day
-
-### Phase 3: Add Derived Fields & Statistics (Important) ✅ COMPLETED
-
-**Priority:** P1 - Useful for post-processing
-
-**Status:** Completed on 2026-01-01
-
-**Tasks:**
-
-- [x] **3.1 Implement field statistics function**
-  - `calculate_field_stats(data)` - Compute min, max, avg, sum for a field
-  - Returns dict with 'min', 'max', 'avg', 'sum' keys
-
-- [x] **3.2 Implement velocity magnitude computation**
-  - `compute_velocity_magnitude(u, v, nx, ny)` - Compute sqrt(u^2 + v^2)
-  - Returns list of velocity magnitudes
-
-- [x] **3.3 Implement comprehensive flow statistics**
-  - `compute_flow_statistics(u, v, p, nx, ny)` - Statistics for all flow components
-  - Returns dict with 'u', 'v', 'p', 'velocity_magnitude' stats
-
-- [x] **3.4 Add tests**
-  - Created `tests/test_derived_fields.py` with comprehensive tests
-  - Tests for all three functions with edge cases
-  - Proper error handling tests (empty lists, wrong types, size mismatches)
-
-**Actual effort:** < 1 day
-
-### Phase 4: Add Error Handling API (Important) ✅ COMPLETED
-
-**Priority:** P1 - Better debugging
-
-**Status:** Completed on 2026-01-02
-
-**Tasks:**
-
-- [x] **4.1 Expose error functions**
-  - `get_last_error()` → Python string (already in C extension)
-  - `get_last_status()` → Python enum (already in C extension)
-  - `get_error_string(code)` → Python string (already in C extension)
-  - `clear_error()` (already in C extension)
-
-- [x] **4.2 Create Python exceptions**
-  - Created `cfd_python/_exceptions.py` with exception hierarchy
-  - `CFDError` base exception with `status_code` and `message` attributes
-  - `CFDMemoryError(CFDError, MemoryError)` - for CFD_ERROR_NOMEM (-2)
-  - `CFDInvalidError(CFDError, ValueError)` - for CFD_ERROR_INVALID (-3)
-  - `CFDIOError(CFDError, IOError)` - for CFD_ERROR_IO (-4)
-  - `CFDUnsupportedError(CFDError, NotImplementedError)` - for CFD_ERROR_UNSUPPORTED (-5)
-  - `CFDDivergedError(CFDError)` - for CFD_ERROR_DIVERGED (-6)
-  - `CFDMaxIterError(CFDError)` - for CFD_ERROR_MAX_ITER (-7)
-
-- [x] **4.3 Implement raise_for_status helper**
-  - `raise_for_status(status_code, context="")` - Raises appropriate exception based on status code
-  - Maps status codes to exception classes
-  - Includes error message from C library when available
-
-- [x] **4.4 Add tests**
-  - Added tests to `tests/test_errors.py`
-  - Tests for exception class hierarchy and inheritance
-  - Tests for `raise_for_status` function with all error codes
-  - Tests for export verification
-
-**Actual effort:** < 0.5 days
-
-### Phase 5: Add Backend Availability API (v0.1.6 Feature) ✅ COMPLETED
-
-**Priority:** P1 - Important for v0.1.6 compatibility
-
-**Status:** Completed on 2025-12-31
-
-**Tasks:**
-
-- [x] **5.1 Expose backend enum**
-  - Added `BACKEND_SCALAR`, `BACKEND_SIMD`, `BACKEND_OMP`, `BACKEND_CUDA` constants
-  - Map to `ns_solver_backend_t` enum (values 0-3)
-
-- [x] **5.2 Implement backend availability functions**
-  - `backend_is_available(backend)` → bool
-  - `backend_get_name(backend)` → string
-  - `list_solvers_by_backend(backend)` → list of solver names
-
-- [x] **5.3 Add backend query helpers**
-  - `get_available_backends()` → list of available backend names
-
-- [x] **5.4 Add tests**
-  - Created `tests/test_backend_availability.py` with comprehensive tests
-  - Tests for constants, availability checking, name queries, solver listing
-
-**Note:** `create_solver_checked()` and `get_solver_backend()` deferred to Phase 4 (error handling integration).
-
-**Actual effort:** 0.5 days
-
-### Phase 6: Add CPU Features & Misc (Enhancement) ✅ COMPLETED
-
-**Priority:** P2 - Nice to have
-
-**Status:** Completed on 2026-01-02
-
-**Tasks:**
-
-- [x] **6.1 CPU features detection**
-  - Added `SIMD_NONE`, `SIMD_AVX2`, `SIMD_NEON` constants
-  - `get_simd_arch()` → returns SIMD architecture constant
-  - `get_simd_name()` → returns "avx2", "neon", or "none"
-  - `has_avx2()` → bool (checks AVX2 availability with OS support verification)
-  - `has_neon()` → bool (checks ARM NEON availability)
-  - `has_simd()` → bool (checks if any SIMD is available)
-
-- [x] **6.2 Grid initialization variants**
-  - `create_grid_stretched(nx, ny, xmin, xmax, ymin, ymax, beta)` - Hyperbolic cosine stretching
-  - Clusters points toward the domain center (symmetric stretching)
-  - Note: Chebyshev and geometric variants not available in C library
-
-- [x] **6.3 Add tests**
-  - Created `tests/test_cpu_features.py` with 26 tests
-  - Tests for SIMD constants, detection functions, and grid stretching
-
-**Actual effort:** < 0.5 days
-
-### Phase 7: Documentation & Tests (Required) ✅ COMPLETED
-
-**Priority:** P1 - Required for release
-
-**Status:** Completed on 2026-01-03
-
-**Tasks:**
-
-- [x] **7.1 Update README**
-  - Complete API reference for all v0.1.6 features
-  - Boundary conditions documentation with examples
-  - Backend availability API documentation
-  - Error handling with exception classes
-  - CPU features detection
-  - Migration guide from v0.1.0
-
-- [x] **7.2 Update Python docstrings**
-  - Comprehensive module docstring in `__init__.py`
-  - All BC functions, types, and edges documented
-  - Error handling API documented
-  - Backend availability API documented
-  - CPU features API documented
-
-- [x] **7.3 Add comprehensive tests**
-  - `test_boundary_conditions.py` - All BC types and backends
-  - `test_backend_availability.py` - Backend constants and functions
-  - `test_derived_fields.py` - Statistics and velocity magnitude
-  - `test_errors.py` - Exception classes and raise_for_status
-  - `test_cpu_features.py` - SIMD detection and grid stretching
-  - `test_abi_compatibility.py` - NULL handling and stress tests
-
-- [x] **7.4 Update examples**
-  - `boundary_conditions.py` - BC types, backends, inlet/outlet
-  - `backend_detection.py` - CPU features, solver/BC backends
-  - `derived_fields.py` - Statistics, velocity magnitude, VTK output
-  - `error_handling.py` - Exception classes, raise_for_status
-
-**Actual effort:** < 0.5 days
-
----
-
-## File Changes Summary
-
-### Files to Modify
-
-| File | Changes |
-|------|---------|
-| `src/cfd_python.c` | Major rewrite - types, registry, errors, new functions |
-| `CMakeLists.txt` | CFD library detection, version check |
-| `cfd_python/__init__.py` | Export new functions, enums, exceptions |
-| `cfd_python/_loader.py` | No changes expected |
-| `pyproject.toml` | Version bump, dependency update |
-| `README.md` | Update documentation |
-
-### Files to Remove
-
-| File | Reason |
-|------|--------|
-| `src/cfd_lib/include/*.h` | Bundled headers - use installed library instead |
-
-### Files to Create
-
-| File | Purpose |
-|------|---------|
-| `src/boundary_conditions.c` | BC wrapper functions (optional - can be in main file) |
-| `tests/test_boundary_conditions.py` | BC tests |
-| `tests/test_derived_fields.py` | Derived fields tests |
-| `tests/test_error_handling.py` | Error handling tests |
-| `tests/test_backend_availability.py` | Backend availability API tests |
-| `examples/boundary_conditions.py` | BC usage examples |
-| `examples/backend_detection.py` | Backend detection examples |
-
----
-
-## Dependency Changes
-
-### Build Dependencies
-
-```toml
-# pyproject.toml
-[build-system]
-requires = ["scikit-build-core", "numpy"]
-
-[project]
-dependencies = ["numpy>=1.20"]
-```
-
-### Runtime Dependencies
-
-- CFD library >= 0.1.6 installed system-wide
-- NumPy >= 1.20
-
-### CMake Requirements
-
-```cmake
-# Find CFD library (v0.1.6 with modular backend libraries)
-find_package(CFD 0.1.6 REQUIRED)
-
-# Link against the unified library (includes all backends)
-target_link_libraries(cfd_python PRIVATE CFD::Library)
-
-# Or manual detection
-find_path(CFD_INCLUDE_DIR cfd/api/simulation_api.h)
-find_library(CFD_LIBRARY cfd_library)  # Unified library name
-```
-
----
-
-## Timeline Estimate
-
-| Phase | Duration | Cumulative |
-|-------|----------|------------|
-| Phase 1: Breaking Changes | ~~2-3 days~~ ✅ 1 day | ~~2-3 days~~ 1 day |
-| Phase 2: Boundary Conditions | ~~3-4 days~~ ✅ 1 day | ~~5-7 days~~ 2 days |
-| Phase 2.5: CI/Build System (v0.1.6) | ✅ 1 day | 3 days |
-| Phase 3: Derived Fields | ~~1-2 days~~ ✅ < 1 day | 3.5 days |
-| Phase 4: Error Handling | ~~1 day~~ ✅ < 0.5 days | 4 days |
-| Phase 5: Backend Availability (v0.1.6) | ✅ 0.5 days | 4.5 days |
-| Phase 6: CPU Features | ~~1 day~~ ✅ < 0.5 days | 5 days |
-| Phase 7: Docs & Tests | ~~2 days~~ ✅ < 0.5 days | 5.5 days |
-
-**Total estimated effort:** ~~9-10 days~~ ~5.5 days ✅ ALL PHASES COMPLETE
-
----
-
-## Risk Assessment
-
-| Risk | Impact | Mitigation |
-|------|--------|------------|
-| API changes in CFD library | High | Pin to specific version, add version checks |
-| Build system complexity | Medium | Test on all platforms in CI |
-| BC backend compatibility | Medium | Test with fallback to scalar backend |
-| NumPy ABI compatibility | Low | Use stable ABI, test multiple versions |
-
----
-
-## Success Criteria
-
-1. All existing tests pass
-2. New BC tests pass
-3. Backend availability API tests pass
-4. Builds successfully on Windows, Linux, macOS
-5. Works with CFD library v0.1.6
-6. Correctly detects available backends at runtime
-7. Python API is Pythonic and well-documented
-8. Examples run successfully
-
----
-
-## Pythonic API Assessment
-
-**Assessment Date:** 2026-01-02
-
-### ✅ What's Pythonic (Good Practices)
-
-**1. Exception Hierarchy Design** - Excellent
-- Multiple inheritance with standard Python exceptions (`CFDMemoryError(CFDError, MemoryError)`)
-- Enables idiomatic exception handling: `except MemoryError` catches CFD memory errors
-- The `raise_for_status()` pattern follows requests library conventions
-- Docstrings with attribute documentation
-
-**2. Naming Conventions** - Good
-- Functions use `snake_case` (`create_grid`, `run_simulation`, `has_avx2`)
-- Constants use `UPPER_CASE` (`SIMD_AVX2`, `BC_TYPE_DIRICHLET`)
-- Private modules prefixed with underscore (`_loader.py`, `_exceptions.py`)
-
-**3. Return Types** - Good
-- Functions return native Python types (dict, list, bool, int) instead of custom objects
-- `create_grid()` returns a dict with intuitive keys (`nx`, `ny`, `x_coords`)
-- `run_simulation_with_params()` returns structured dict with `stats` sub-dict
-
-**4. Module Organization** - Good
-- Clean `__all__` exports
-- Graceful handling of unbuilt extension (`ExtensionNotBuiltError`)
-- Dynamic solver constants discovery
-
-### ⚠️ Areas for Improvement
-
-**1. ~~Inconsistent Coordinate Naming~~ (FIXED)**
-```python
-# Both functions now return consistent keys:
-grid["x_coords"], grid["y_coords"]  # standardized naming
-```
-
-**2. Mixed Paradigms for Complex Operations**
-```python
-# Positional args (C-style):
-create_grid(10, 10, 0.0, 1.0, 0.0, 1.0)
-
-# More Pythonic alternative:
-create_grid(nx=10, ny=10, bounds=(0.0, 1.0, 0.0, 1.0))
-# or
-create_grid(nx=10, ny=10, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0)
-```
-
-**3. Verbose Boundary Condition API**
-The BC functions require separate calls for each boundary type. A more Pythonic approach might be:
-```python
-# Current:
-bc_apply_inlet_uniform(u, v, nx, ny, u_inlet, v_inlet, BC_EDGE_LEFT)
-bc_apply_outlet_velocity(u, v, nx, ny, BC_EDGE_RIGHT)
-
-# More Pythonic (context manager or builder pattern):
-with BoundaryConditions(u, v, nx, ny) as bc:
-    bc.left.inlet_uniform(u=1.0, v=0.0)
-    bc.right.outlet()
-```
-
-**4. Constants Could Use IntEnum**
-```python
-# Current:
-SIMD_NONE = 0
-SIMD_AVX2 = 1
-
-# More Pythonic:
-class SIMDArch(IntEnum):
-    NONE = 0
-    AVX2 = 1
-    NEON = 2
-```
-This provides `str(SIMDArch.AVX2)` → `"SIMDArch.AVX2"` and better IDE support.
-
-**5. No Type Hints in Public API**
-While `_exceptions.py` has type hints, the C extension functions lack type stubs (`.pyi` files) for IDE autocompletion.
-
-### 📊 Summary
-
-| Aspect | Rating | Notes |
-|--------|--------|-------|
-| Naming | ★★★★☆ | Follows PEP 8, minor inconsistencies |
-| Error Handling | ★★★★★ | Excellent exception hierarchy |
-| Return Types | ★★★★☆ | Native types, dict structure could use dataclasses |
-| API Design | ★★★☆☆ | Functional but verbose, mirrors C API closely |
-| Type Annotations | ★★☆☆☆ | Missing .pyi stubs for C extension |
-| Documentation | ★★★★☆ | Good docstrings in `__init__.py` |
-
-**Overall:** The code is reasonably Pythonic for a C extension binding. The exception handling is particularly well-designed. The main improvements would be adding type stubs and potentially providing higher-level Pythonic wrappers around the lower-level C-style functions.
-
-### Recommended Future Enhancements
-
-1. **Add type stubs** (`cfd_python.pyi`) for IDE autocompletion
-2. **Standardize coordinate naming** across grid functions
-3. **Consider IntEnum** for constant groups (SIMD, BC types, backends)
-4. **Optional high-level wrappers** for BC operations (Phase 8 candidate)
-
----
-
-## Phase 7 (CFD ROADMAP) Status: Python Integration
-
-This section tracks progress on Phase 7 from the main CFD library ROADMAP.md, which covers Python integration.
-
-### 7.1 Python Bindings (P1)
-
-| Task | Status | Notes |
-|------|--------|-------|
-| Complete C extension module | ✅ Done | `cfd_python.c` with Stable ABI (Py_LIMITED_API 3.9+) |
-| NumPy array integration | ⚠️ Partial | Functions accept/return Python lists; NumPy conversion at Python layer |
-| Pythonic API design | ✅ Done | Snake_case naming, native Python types, exception hierarchy |
-| Type hints and stubs | ❌ Pending | Phase 8 in cfd-python ROADMAP - `.pyi` stubs not yet created |
-| Pre-built wheels (manylinux, macOS, Windows) | ✅ Done | CI builds wheels for all platforms; CUDA variant for Linux/Windows |
-
-**C Extension Module Features:**
-- Grid creation (uniform and stretched)
-- Flow field operations
-- Solver registry and creation
-- Simulation API (`run_simulation`, `run_simulation_step`)
-- Boundary conditions (all types: Neumann, Dirichlet, no-slip, inlet, outlet)
-- Derived fields and statistics
-- Error handling with Python exceptions
-- Backend availability API (Scalar, SIMD, OMP, CUDA)
-- CPU features detection (AVX2, NEON)
-
-### 7.2 High-level Python API (P1)
-
-| Task | Status | Notes |
-|------|--------|-------|
-| Problem definition classes | ❌ Pending | Not yet implemented |
-| Mesh generation helpers | ⚠️ Partial | `create_grid`, `create_grid_stretched` available |
-| Post-processing utilities | ⚠️ Partial | `compute_velocity_magnitude`, `compute_flow_statistics` available |
-| Jupyter notebook integration | ❌ Pending | Not yet implemented |
-
-**What's Available:**
-- Low-level C bindings exposed directly to Python
-- Exception hierarchy for error handling (`CFDError` and subclasses)
-- `raise_for_status()` helper for checking C library return codes
-
-**What's Missing:**
-- High-level `Simulation` class for problem setup
-- `BoundaryConditions` builder/context manager pattern
-- Integration with matplotlib/visualization in notebooks
-- Problem templates (lid-driven cavity, channel flow, etc.)
-
-### Summary
-
-Phase 7.1 (Python Bindings) is **largely complete** - the C extension module is functional, builds wheels for all platforms, and provides a working Python API. The main gaps are:
-1. Type stubs for IDE support (Phase 8 in cfd-python)
-2. NumPy direct integration (currently uses list conversion)
-
-Phase 7.2 (High-level Python API) is **not started** - the current API exposes C functions directly without high-level Pythonic wrappers or classes.
diff --git a/ROADMAP.md b/ROADMAP.md
index 3a50e5a..6a2d0f9 100644
--- a/ROADMAP.md
+++ b/ROADMAP.md
@@ -2,8 +2,8 @@
 
 This document outlines planned enhancements and new features for cfd-python after the v0.1.6 migration is complete.
 
-**Last Updated:** 2026-01-02
-**Current Version:** 0.1.6 (migration in progress)
+**Last Updated:** 2026-03-05
+**Current Version:** 0.1.6 (released)
 **Target Version:** 0.2.0+
 
 ---
@@ -14,10 +14,70 @@ With the v0.1.6 migration completing core functionality (boundary conditions, er
 
 ---
 
-## Phase 8: Type Stubs & IDE Support
+## CI/CD Lessons Learned (v0.1.6 Release)
+
+Key insights from the PyPI publishing process that should inform future releases:
+
+### Manylinux Compatibility (Critical for PyPI)
+
+| Issue | Solution |
+|-------|----------|
+| PyPI rejects `linux_x86_64` wheels | Build inside manylinux containers |
+| glibc version mismatch on ubuntu-latest | Use `quay.io/pypa/manylinux_2_28_x86_64` container |
+| auditwheel requires patchelf >= 0.14 | Install via pip (EPEL version too old) |
+| CUDA wheels need manylinux too | Use `nvidia/cuda:12.4.0-devel-rockylinux8` (manylinux_2_28 compatible) |
+
+### CUDA Wheel Building
+
+| Issue | Solution |
+|-------|----------|
+| CUDA not in standard manylinux containers | Use NVIDIA CUDA devel images based on Rocky Linux 8 |
+| CMake FetchContent needs git | `dnf install -y git` in container |
+| patchelf too old in EPEL | `pip install patchelf` instead |
+
+### Version Management
+
+| Issue | Solution |
+|-------|----------|
+| PyPI rejects local versions (`+g...`) | Build from git tag, not branch |
+| workflow_dispatch defaults to branch | Add `ref` input parameter to workflow |
+| sdist filename with hyphens | PyPI requires underscores (PEP 625) - rename if needed |
+
+### Partial Upload Recovery
+
+| Issue | Solution |
+|-------|----------|
+| Upload fails midway, can't re-upload same files | Add `skip-existing: true` to pypi-publish action |
+
+### Windows CUDA Optimization
+
+| Issue | Solution |
+|-------|----------|
+| CUDA toolkit install takes 14-16 min | Use `method: network` with minimal sub-packages |
+| CMake can't find CUDA toolset | Include `visual_studio_integration` in sub-packages |
+
+### Architecture Support
+
+For future ARM64 Linux support:
+
+```bash
+ARCH=$(uname -m)
+auditwheel repair dist_raw/*.whl --plat manylinux_2_28_${ARCH} -w dist/
+```
+
+### Recommended Workflow Structure
+
+1. **Build C library and wheel inside manylinux container** (not on host)
+2. **Use auditwheel to repair and tag wheels**
+3. **Support tag-based builds via `ref` input for releases**
+4. **Always use `skip-existing: true` for PyPI uploads**
+
+---
+
+## Phase 8: Type Stubs & IDE Support ✓ (Mostly Complete)
 
 **Priority:** P1 - High impact for developer experience
-**Estimated Effort:** 2-3 days
+**Status:** Core tasks done, CI integration remaining
 
 ### Goals
 
@@ -27,60 +87,23 @@ With the v0.1.6 migration completing core functionality (boundary conditions, er
 
 ### Tasks
 
-- [ ] **8.1 Create `cfd_python.pyi` stub file**
-  ```python
-  # Example stub content
-  from typing import Dict, List, Optional, Union
-
-  # Constants
-  SIMD_NONE: int
-  SIMD_AVX2: int
-  SIMD_NEON: int
-
-  BACKEND_SCALAR: int
-  BACKEND_SIMD: int
-  BACKEND_OMP: int
-  BACKEND_CUDA: int
-
-  # Functions
-  def create_grid(
-      nx: int,
-      ny: int,
-      xmin: float,
-      xmax: float,
-      ymin: float,
-      ymax: float,
-  ) -> Dict[str, Union[int, float, List[float]]]: ...
-
-  def run_simulation(
-      nx: int,
-      ny: int,
-      *,
-      steps: int = ...,
-      solver_type: str = ...,
-      xmin: float = ...,
-      xmax: float = ...,
-      ymin: float = ...,
-      ymax: float = ...,
-      output_file: Optional[str] = ...,
-  ) -> List[float]: ...
-  ```
+- [x] **8.1 Create `cfd_python/__init__.pyi` stub file**
+  - Comprehensive stubs covering all exported functions, constants, and classes
 
-- [ ] **8.2 Add `py.typed` marker file**
-  - Create empty `cfd_python/py.typed` for PEP 561 compliance
+- [x] **8.2 Add `py.typed` marker file**
+  - `cfd_python/py.typed` exists for PEP 561 compliance
 
-- [ ] **8.3 Update `pyproject.toml`**
-  - Add stub files to package data
-  - Configure mypy/pyright settings
+- [x] **8.3 Update `pyproject.toml`**
+  - mypy included in `[project.optional-dependencies] dev`
 
 - [ ] **8.4 Add type checking to CI**
   - Run mypy on tests to verify stubs are correct
 
 ### Success Criteria
 
-- IDE shows function signatures and return types
-- `mypy tests/` passes without errors
-- Stubs cover all exported functions and constants
+- [x] IDE shows function signatures and return types
+- [ ] `mypy tests/` passes without errors — not yet in CI
+- [x] Stubs cover all exported functions and constants
 
 ---
 
@@ -98,66 +121,6 @@ With the v0.1.6 migration completing core functionality (boundary conditions, er
 ### Tasks
 
 - [ ] **9.1 Create enum classes in `_enums.py`**
-  ```python
-  from enum import IntEnum
-
-  class SIMDArch(IntEnum):
-      """SIMD architecture detected on the system."""
-      NONE = 0
-      AVX2 = 1
-      NEON = 2
-
-  class Backend(IntEnum):
-      """Solver backend types."""
-      SCALAR = 0
-      SIMD = 1
-      OMP = 2
-      CUDA = 3
-
-  class BCType(IntEnum):
-      """Boundary condition types."""
-      PERIODIC = 0
-      NEUMANN = 1
-      DIRICHLET = 2
-      NOSLIP = 3
-      INLET = 4
-      OUTLET = 5
-
-  class BCEdge(IntEnum):
-      """Boundary edges."""
-      LEFT = 0
-      RIGHT = 1
-      BOTTOM = 2
-      TOP = 3
-
-  class BCBackend(IntEnum):
-      """Boundary condition backends."""
-      AUTO = 0
-      SCALAR = 1
-      OMP = 2
-      SIMD = 3
-      CUDA = 4
-
-  class OutputType(IntEnum):
-      """Output file types."""
-      VELOCITY = 0
-      VELOCITY_MAGNITUDE = 1
-      FULL_FIELD = 2
-      CSV_TIMESERIES = 3
-      CSV_CENTERLINE = 4
-      CSV_STATISTICS = 5
-
-  class StatusCode(IntEnum):
-      """CFD operation status codes."""
-      SUCCESS = 0
-      ERROR = -1
-      ERROR_NOMEM = -2
-      ERROR_INVALID = -3
-      ERROR_IO = -4
-      ERROR_UNSUPPORTED = -5
-      ERROR_DIVERGED = -6
-      ERROR_MAX_ITER = -7
-  ```
 
 - [ ] **9.2 Export enums alongside bare constants**
   - Maintain backward compatibility with `SIMD_NONE`, `BACKEND_SCALAR`, etc.
@@ -191,169 +154,12 @@ With the v0.1.6 migration completing core functionality (boundary conditions, er
 ### Tasks
 
 - [ ] **10.1 Create `Grid` class**
-  ```python
-  class Grid:
-      """High-level grid object with coordinate access."""
-
-      def __init__(self, nx: int, ny: int,
-                   xmin: float = 0.0, xmax: float = 1.0,
-                   ymin: float = 0.0, ymax: float = 1.0,
-                   stretching: Optional[float] = None):
-          if stretching is not None:
-              self._data = create_grid_stretched(nx, ny, xmin, xmax, ymin, ymax, stretching)
-          else:
-              self._data = create_grid(nx, ny, xmin, xmax, ymin, ymax)
-
-      @property
-      def nx(self) -> int:
-          return self._data["nx"]
-
-      @property
-      def ny(self) -> int:
-          return self._data["ny"]
-
-      @property
-      def x(self) -> List[float]:
-          return self._data.get("x", self._data.get("x_coords", []))
-
-      @property
-      def y(self) -> List[float]:
-          return self._data.get("y", self._data.get("y_coords", []))
-
-      @property
-      def shape(self) -> tuple[int, int]:
-          return (self.nx, self.ny)
-
-      @property
-      def bounds(self) -> tuple[float, float, float, float]:
-          return (self._data["xmin"], self._data["xmax"],
-                  self._data["ymin"], self._data["ymax"])
-  ```
 
 - [ ] **10.2 Create `BoundaryConditions` builder**
-  ```python
-  class BoundaryConditions:
-      """Fluent builder for boundary conditions."""
-
-      def __init__(self, u: List[float], v: List[float], nx: int, ny: int):
-          self.u = u
-          self.v = v
-          self.nx = nx
-          self.ny = ny
-
-      def noslip_walls(self) -> "BoundaryConditions":
-          """Apply no-slip condition to all walls."""
-          bc_apply_noslip(self.u, self.v, self.nx, self.ny)
-          return self
-
-      def inlet_uniform(self, edge: BCEdge, u: float, v: float = 0.0) -> "BoundaryConditions":
-          """Apply uniform inlet velocity."""
-          bc_apply_inlet_uniform(self.u, self.v, self.nx, self.ny, u, v, int(edge))
-          return self
-
-      def inlet_parabolic(self, edge: BCEdge, max_velocity: float) -> "BoundaryConditions":
-          """Apply parabolic inlet profile."""
-          bc_apply_inlet_parabolic(self.u, self.v, self.nx, self.ny, max_velocity, int(edge))
-          return self
-
-      def outlet(self, edge: BCEdge) -> "BoundaryConditions":
-          """Apply zero-gradient outlet."""
-          bc_apply_outlet_velocity(self.u, self.v, self.nx, self.ny, int(edge))
-          return self
-
-      # Context manager support
-      def __enter__(self) -> "BoundaryConditions":
-          return self
-
-      def __exit__(self, *args) -> None:
-          pass
-
-  # Usage:
-  # bc = BoundaryConditions(u, v, nx, ny)
-  # bc.noslip_walls().inlet_uniform(BCEdge.LEFT, u=1.0).outlet(BCEdge.RIGHT)
-  ```
 
 - [ ] **10.3 Create `Simulation` class**
-  ```python
-  class Simulation:
-      """High-level simulation runner."""
-
-      def __init__(self, grid: Grid, solver: str = "explicit_euler"):
-          self.grid = grid
-          self.solver = solver
-          self._params = get_default_solver_params()
-          self._result = None
-
-      def set_params(self, **kwargs) -> "Simulation":
-          """Set solver parameters."""
-          self._params.update(kwargs)
-          return self
-
-      def run(self, steps: int, output_file: Optional[str] = None) -> "SimulationResult":
-          """Run the simulation."""
-          result = run_simulation_with_params(
-              self.grid.nx, self.grid.ny,
-              *self.grid.bounds,
-              steps=steps,
-              solver_type=self.solver,
-              output_file=output_file,
-              **self._params
-          )
-          self._result = SimulationResult(result)
-          return self._result
-
-  class SimulationResult:
-      """Container for simulation results with statistics."""
-
-      def __init__(self, data: dict):
-          self._data = data
-
-      @property
-      def velocity_magnitude(self) -> List[float]:
-          return self._data["velocity_magnitude"]
-
-      @property
-      def stats(self) -> dict:
-          return self._data["stats"]
-
-      def compute_statistics(self, u: List[float], v: List[float], p: List[float]) -> dict:
-          return compute_flow_statistics(u, v, p, self._data["nx"], self._data["ny"])
-  ```
 
 - [ ] **10.4 Create `Field` class for data manipulation**
-  ```python
-  class Field:
-      """Scalar field with statistics and operations."""
-
-      def __init__(self, data: List[float], nx: int, ny: int):
-          self.data = data
-          self.nx = nx
-          self.ny = ny
-
-      @property
-      def shape(self) -> tuple[int, int]:
-          return (self.nx, self.ny)
-
-      @property
-      def stats(self) -> dict:
-          return calculate_field_stats(self.data)
-
-      @property
-      def min(self) -> float:
-          return self.stats["min"]
-
-      @property
-      def max(self) -> float:
-          return self.stats["max"]
-
-      @property
-      def mean(self) -> float:
-          return self.stats["avg"]
-
-      def to_vtk(self, filename: str, name: str = "field") -> None:
-          write_vtk_scalar(filename, name, self.data, self.nx, self.ny,
-                          0.0, 1.0, 0.0, 1.0)
-  ```
 
 - [ ] **10.5 Add tests for high-level API**
 
@@ -385,6 +191,7 @@ With the v0.1.6 migration completing core functionality (boundary conditions, er
   - Support both list and ndarray inputs
 
 - [ ] **11.2 Create array conversion utilities**
+
   ```python
   def to_numpy(data: List[float], shape: tuple[int, int]) -> np.ndarray:
       """Convert flat list to 2D NumPy array."""
@@ -396,6 +203,7 @@ With the v0.1.6 migration completing core functionality (boundary conditions, er
   ```
 
 - [ ] **11.3 Add `as_numpy` option to functions**
+
   ```python
   def run_simulation(..., as_numpy: bool = False) -> Union[List[float], np.ndarray]:
       result = _run_simulation_impl(...)
@@ -405,6 +213,7 @@ With the v0.1.6 migration completing core functionality (boundary conditions, er
   ```
 
 - [ ] **11.4 Support array protocol in Field class**
+
   ```python
   class Field:
       def __array__(self) -> np.ndarray:
@@ -439,6 +248,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
 ### Tasks
 
 - [ ] **12.1 Add optional dependency**
+
   ```toml
   [project.optional-dependencies]
   viz = ["cfd-visualization>=0.2.0"]
@@ -469,6 +279,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
 ### Tasks
 
 - [ ] **13.1 Add timing to simulation results**
+
   ```python
   result = run_simulation_with_params(...)
   print(result["timing"])
@@ -476,6 +287,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
   ```
 
 - [ ] **13.2 Create benchmarking utilities**
+
   ```python
   from cfd_python.benchmark import benchmark_solver
 
@@ -489,6 +301,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
   ```
 
 - [ ] **13.3 Add backend comparison utility**
+
   ```python
   from cfd_python.benchmark import compare_backends
 
@@ -511,6 +324,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
 
 **Priority:** P3 - Advanced use case
 **Estimated Effort:** 3-5 days
+**Spec:** [`.claude/specs/phase-14-async-parallel-simulation.md`](.claude/specs/phase-14-async-parallel-simulation.md)
 
 ### Goals
 
@@ -520,39 +334,19 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
 
 ### Tasks
 
-- [ ] **14.1 Add progress callback support**
-  ```python
-  def progress_callback(step: int, total: int, stats: dict):
-      print(f"Step {step}/{total}: max_vel={stats['max_velocity']:.4f}")
-
-  run_simulation_with_params(..., callback=progress_callback, callback_interval=10)
-  ```
-
-- [ ] **14.2 Create async simulation wrapper**
-  ```python
-  import asyncio
-
-  async def run_simulation_async(...) -> SimulationResult:
-      """Run simulation in background thread."""
-      loop = asyncio.get_event_loop()
-      return await loop.run_in_executor(None, run_simulation_with_params, ...)
-  ```
-
-- [ ] **14.3 Add parameter sweep utility**
-  ```python
-  from cfd_python.parallel import parameter_sweep
-
-  results = parameter_sweep(
-      base_params={'nx': 64, 'ny': 64, 'steps': 100},
-      vary={'dt': [0.001, 0.0005, 0.0001], 'cfl': [0.5, 0.8]},
-      n_workers=4
-  )
-  ```
+- [ ] **14.1 Add progress callback support** — `callback` and `callback_interval` keyword arguments on `run_simulation_with_params()`, with cancellation via `False` return
+- [ ] **14.1b Release GIL during simulation loop** — `Py_BEGIN_ALLOW_THREADS` / `Py_END_ALLOW_THREADS` around the step loop to enable thread-based parallelism
+- [ ] **14.2 Create async simulation wrapper** — `run_simulation_async()` in `cfd_python/_async.py` using `loop.run_in_executor()` with task cancellation support
+- [ ] **14.3 Add parameter sweep utility** — `parameter_sweep()` in `cfd_python/_parallel.py` using `ProcessPoolExecutor` with Cartesian product of parameter variations
+- [ ] **14.4 Update exports and type stubs** — re-export new functions from `__init__.py`, add signatures to `__init__.pyi`
+- [ ] **14.5 Add tests** — `tests/test_async.py` and `tests/test_parallel.py` covering callbacks, cancellation, concurrency, and parameter sweeps
 
 ### Success Criteria
 
-- Long simulations can report progress
-- Multiple simulations can run in parallel
+- Long simulations can report progress via callbacks without significant overhead
+- Callbacks can cancel a running simulation
+- `run_simulation_async()` does not block the asyncio event loop
+- Multiple simulations run in parallel via `parameter_sweep()` with true CPU parallelism
 - Parameter sweeps are easy to set up
 
 ---
@@ -571,6 +365,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
 ### Tasks
 
 - [ ] **15.1 Expose 3D grid functions from C library**
+
   ```python
   def create_grid_3d(
       nx: int, ny: int, nz: int,
@@ -581,6 +376,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
   ```
 
 - [ ] **15.2 Add 3D boundary condition edges**
+
   ```python
   class BCFace(IntEnum):
       """3D boundary faces."""
@@ -593,6 +389,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
   ```
 
 - [ ] **15.3 Extend Field class for 3D**
+
   ```python
   class Field3D(Field):
       def __init__(self, data: List[float], nx: int, ny: int, nz: int):
@@ -630,6 +427,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
 ### Tasks
 
 - [ ] **16.1 Add checkpoint save/load functions**
+
   ```python
   def save_checkpoint(
       filename: str,
@@ -646,6 +444,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
   ```
 
 - [ ] **16.2 Add automatic checkpointing to Simulation class**
+
   ```python
   sim = Simulation(grid)
   sim.run(
@@ -656,12 +455,14 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
   ```
 
 - [ ] **16.3 Add resume functionality**
+
   ```python
   sim = Simulation.from_checkpoint("./checkpoints/step_5000.h5")
   sim.run(steps=5000)  # Continues from step 5000
   ```
 
 - [ ] **16.4 Optional HDF5 support**
+
   ```toml
   [project.optional-dependencies]
   hdf5 = ["h5py>=3.0"]
@@ -689,6 +490,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
 ### Tasks
 
 - [ ] **17.1 Implement lid-driven cavity benchmark**
+
   ```python
   from cfd_python.validation import lid_driven_cavity
 
@@ -701,6 +503,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
   ```
 
 - [ ] **17.2 Implement Poiseuille flow validation**
+
   ```python
   from cfd_python.validation import poiseuille_flow
 
@@ -712,6 +515,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) for visualizati
   ```
 
 - [ ] **17.3 Implement Taylor-Green vortex decay**
+
   ```python
   from cfd_python.validation import taylor_green_vortex
 
@@ -822,6 +626,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
 ### Tasks
 
 - [ ] **19.1 Define configuration schema**
+
   ```yaml
   # simulation.yaml
   grid:
@@ -865,6 +670,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
   ```
 
 - [ ] **19.2 Create config loader**
+
   ```python
   from cfd_python.config import load_config, run_from_config
 
@@ -873,10 +679,11 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
   ```
 
 - [ ] **19.3 Add CLI entry point**
+
   ```bash
-  $ cfd-python run simulation.yaml
-  $ cfd-python validate simulation.yaml
-  $ cfd-python info  # Show available solvers, backends
+  cfd-python run simulation.yaml
+  cfd-python validate simulation.yaml
+  cfd-python info  # Show available solvers, backends
   ```
 
 - [ ] **19.4 Config validation with helpful errors**
@@ -903,6 +710,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
 ### Tasks
 
 - [ ] **20.1 Define plugin interface**
+
   ```python
   from cfd_python.plugins import SolverPlugin, register_plugin
 
@@ -917,6 +725,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
   ```
 
 - [ ] **20.2 Plugin discovery mechanism**
+
   ```python
   # Automatic discovery via entry points
   # pyproject.toml of plugin package:
@@ -925,6 +734,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
   ```
 
 - [ ] **20.3 Custom BC plugin support**
+
   ```python
   from cfd_python.plugins import BCPlugin
 
@@ -960,6 +770,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
 ### Tasks
 
 - [ ] **21.1 Add memory usage reporting**
+
   ```python
   from cfd_python.memory import estimate_memory, get_memory_usage
 
@@ -972,6 +783,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
   ```
 
 - [ ] **21.2 Memory-mapped field storage**
+
   ```python
   grid = Grid(nx=4096, ny=4096)
   sim = Simulation(grid, memory_mode="mmap", mmap_dir="/tmp/cfd")
@@ -1005,6 +817,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
 ### Tasks
 
 - [ ] **22.1 GPU memory reporting**
+
   ```python
   from cfd_python.cuda import get_gpu_info, get_gpu_memory
 
@@ -1016,6 +829,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
   ```
 
 - [ ] **22.2 Device selection**
+
   ```python
   from cfd_python.cuda import set_device
 
@@ -1024,6 +838,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
   ```
 
 - [ ] **22.3 Multi-GPU domain decomposition**
+
   ```python
   sim = Simulation(
       grid,
@@ -1057,6 +872,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
 ### Tasks
 
 - [ ] **23.1 Add structured logging**
+
   ```python
   import logging
   logging.basicConfig(level=logging.INFO)
@@ -1067,6 +883,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
   ```
 
 - [ ] **23.2 Configurable log levels**
+
   ```python
   import cfd_python
   cfd_python.set_log_level("DEBUG")  # Verbose output
@@ -1074,6 +891,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
   ```
 
 - [ ] **23.3 Metrics collection**
+
   ```python
   from cfd_python.metrics import get_metrics
 
@@ -1149,6 +967,7 @@ See [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) Phase 2 for vis
 ### Background
 
 Physics-Informed Neural Networks embed physical laws (Navier-Stokes equations) into the loss function, achieving:
+
 - ~1000× speedup over traditional CFD for inference
 - <5% relative error for velocity fields
 - 5-10× less memory than CFD solvers
@@ -1277,7 +1096,7 @@ ml-jax = ["jax>=0.4", "flax>=0.7", "cfd-python[ml]"]
 
 | Version | Phases | Focus |
 |---------|--------|-------|
-| 0.2.0 | 8, 9 | Type safety & IDE support |
+| 0.2.0 | 8 (mostly done), 9 | Type safety & IDE support |
 | 0.3.0 | 10, 11 | High-level API & NumPy |
 | 0.4.0 | 12, 13 | cfd-visualization integration & profiling |
 | 0.5.0 | 14, 17 | Async, parallel & validation |
@@ -1294,6 +1113,7 @@ ml-jax = ["jax>=0.4", "flax>=0.7", "cfd-python[ml]"]
 Items not yet planned but worth considering:
 
 ### Simulation Features
+
 - [ ] Turbulence models (k-ε, k-ω, LES)
 - [ ] Multiphase flow support
 - [ ] Compressible flow solvers
@@ -1302,6 +1122,7 @@ Items not yet planned but worth considering:
 - [ ] Moving mesh support
 
 ### I/O & Interoperability
+
 - [ ] OpenFOAM mesh import/export
 - [ ] CGNS format support
 - [ ] ParaView Catalyst integration
@@ -1309,6 +1130,7 @@ Items not yet planned but worth considering:
 - [ ] STL geometry import
 
 ### Performance
+
 - [ ] Mixed precision (FP32/FP64)
 - [ ] Sparse matrix solvers
 - [ ] Multigrid preconditioners
@@ -1340,7 +1162,6 @@ See [CONTRIBUTING.md](CONTRIBUTING.md) for guidelines.
 
 ## Related Documents
 
-- [MIGRATION_PLAN.md](MIGRATION_PLAN.md) - Current v0.1.6 migration status
 - [README.md](README.md) - User documentation
 - [CHANGELOG.md](CHANGELOG.md) - Version history
 - [cfd-visualization ROADMAP](../cfd-visualization/ROADMAP.md) - Visualization library roadmap