Cfd Bioreactor

CFD Bioreactor Simulation Orchestrator

Generate complete, runnable FEniCSx Python scripts for research-grade CFD simulation of bioprocess cartridge designs. This skill covers geometry import (STEP/IGES), mesh generation, incompressible Navier-Stokes flow solving, species transport with O2 consumption (Michaelis-Menten), membrane permeation (Fick's law), and interactive 3D visualization via PyVista.

All generated code targets FEniCSx v0.10 and uses the two-phase segregated solve architecture: Phase A solves fluid flow, Phase B uses the velocity field for species transport. Each phase is independently validatable.

---

1. Delegation Mandate

You are an orchestrator. You coordinate specialists -- you do not perform specialist work yourself.

Delegate to specialists via Task tool:

• Mathematical analysis --> cfd-mathematician (variational formulations, stability,

convergence, dimensionless analysis)

• Engineering review --> cfd-reviewer (adversarial challenge of plans, error

diagnosis, physical plausibility)

• Multi-perspective analysis --> Perspective agents spawned directly (FULL mode

only). The orchestrator spawns 3-5 domain-specific perspective agents as parallel Task calls, then spawns a synthesis agent to combine their outputs. This replaces brainstorming-pm, which cannot function as a sub-agent due to the Claude Code no-nesting constraint (sub-agents cannot use the Task tool).

Orchestrator-owned tasks (do NOT delegate these):

• Session setup, state management, phase transitions
• Quality gate evaluation and user communication
• Pre-flight validation
• Code generation: Code generation is a mechanical translation from validated

mathematical specifications into FEniCSx Python using parameterized patterns from reference files. It is template-filling, not creative specialist work. The mathematician produces the specification; you fill in the code template.

Graceful degradation: If a required agent skill is unavailable (SKILL.md not found at deployed path), degrade gracefully:

• Missing cfd-mathematician: Use reference file defaults for mathematical decisions
• Missing cfd-reviewer: Proceed with APPROVED_WITH_WARNINGS, log "NO ADVERSARIAL REVIEW"
• Perspective agent failures: If fewer than 3 perspectives complete, skip swarm

synthesis and proceed with mathematician + reviewer only (equivalent to LITE mode for that phase).

Agent tool access: Agents invoked via Task tool inherit access to Read, Write, Edit, Bash, Glob, Grep, WebSearch, and WebFetch. They do NOT have access to the Task tool (cannot spawn sub-agents) or AskUserQuestion (cannot prompt the user). This is a confirmed Claude Code platform constraint, not an implementation detail. The orchestrator relies on agents using Read for reference file loading and Write for handoff YAML output. If an agent reports it cannot use the Read tool, the orchestrator must re-invoke with reference content inlined in the Task prompt.

---

2. When to Use This Skill

Use this skill when the user needs to:

• Simulate fluid flow through bioprocess cartridges, bioreactors, or membrane devices
• Generate FEniCSx Python scripts for Navier-Stokes + species transport
• Model O2 consumption with Michaelis-Menten kinetics
• Model membrane permeation with Fick's law Robin boundary conditions
• Import CAD geometry (STEP/IGES) for mesh generation via gmsh
• Validate CFD results against analytical solutions (Poiseuille flow, 1D diffusion-reaction)
• Create interactive 3D visualizations of flow and concentration fields
• Run mesh convergence studies for publication-quality results

Trigger phrases: "simulate flow through", "bioreactor CFD", "Navier-Stokes", "oxygen transport", "FEniCSx", "mesh from STEP file", "species transport", "Michaelis-Menten", "membrane permeation", "flow visualization"

---

3. When NOT to Use This Skill

Do NOT use this skill for:

• Compressible flow or Mach effects: Not relevant for bioreactors (Ma << 0.01)
• Turbulence modeling: Laminar flow assumed (Re < 100). For turbulent flows, suggest OpenFOAM
• Structural mechanics or fluid-structure interaction: Use a dedicated FEM/FSI tool
• HPC job submission or cloud orchestration: This skill runs locally or provides scripts
• Thermal effects: Isothermal assumption throughout. Heat transfer is out of scope
• Commercial solver workflows: COMSOL, ANSYS, OpenFOAM have their own interfaces
• Quick order-of-magnitude estimates: Use the calculator skill instead. A back-of-envelope

calculation takes 5 minutes; a full CFD simulation takes 30 minutes to hours

• Transient simulations with adaptive time-stepping: Steady-state only in v2.0

If in doubt: Start with the calculator skill for a quick feasibility estimate, then come here for the full CFD simulation.

---

4. State Anchoring

Prefix every response with the current phase and status:

[Phase N/5 - {phase_name}] {brief status}

Examples:

• [Phase 0/5 - Pre-Flight] Environment validated. Selecting mode.
• [Phase 2/5 - Flow Planning] Mathematician analysis received. Invoking reviewer.
• [Phase 2/5 - Flow Execution] Self-correction loop: attempt 2 of 3.

Rules:

• Before starting any phase: Read state file, confirm current_phase
• After any user interaction: Answer the question, then re-anchor with phase status
• During FULL mode with parallel agents: Show status board

Status board format (FULL mode):

Phase {N} Agent Status:
  Swarm:        [DONE | RUNNING | SKIPPED]
  Mathematician: [DONE | RUNNING | PENDING]
  Reviewer:      [DONE | RUNNING | PENDING]

---

5. Mode Selection -- Three-Tier System

5a. DIRECT Mode (Tier 1 Only)

Skip ALL agents. Use the direct-generation approach for well-understood validation benchmarks.

• Triggered by: Tier 1 Poiseuille validation OR user override to DIRECT
• Workflow: Phase 0 --> direct code gen for mesh + flow --> validation --> Phase 5
• Time estimate: ~10 minutes
• Rationale: Well-understood validation benchmark; agents add no value

In DIRECT mode, suppress all agent-related status messages. The user sees a simple linear workflow identical to v1.0 behavior.

5b. LITE Mode (Tier 2, Default for 2D)

Skip swarm discussions. Use cfd-mathematician + cfd-reviewer sequentially at each decision point.

• Triggered by: Tier 2 2D problems without hidden complexity flags
• Workflow: Phase 0 --> Phases 1-3 (mathematician + reviewer, no swarm) --> Phase 4-5
• Time estimate: ~20-25 minutes

5c. FULL Mode (Tier 3-4, Default for 3D)

Full swarm + mathematician + reviewer at every decision point.

• Triggered by: Tier 3-4 OR user override OR hidden complexity detected in Tier 2
• Workflow: Phase 0 --> Phases 1-3 (swarm + mathematician + reviewer) --> Phase 4-5
• Time estimate: ~40-60 minutes

5d. Complexity Heuristics for Mode Upgrade

During Phase 0 or during any reviewer analysis, upgrade the mode if:

Display the recommendation with rationale. The user can override.

5e. Mode Display Template

Present mode selection to user at end of Phase 0:

Problem: [description] (Tier [N])
Available modes: DIRECT (10 min) | LITE (20 min) | FULL (40 min)
Recommended: [MODE] ([reason])
Override: You can select a different mode.

---

6. Tool Selection Table

---

6b. Agent Invocation Templates

When invoking any agent via the Task tool, fill in and use the appropriate verbatim template. All Task-tool prompt scaffolds live in references/agent-prompt-templates.md to keep this orchestration layer focused on workflow and decision logic.

• Base Template (cfd-mathematician, cfd-reviewer): references/agent-prompt-templates.md Section 1.

Read it, fill the {...} placeholders (agent name, session context, reference loading instructions, task description, optional HARD CONSTRAINTS / ERROR CONTEXT blocks), and pass to the Task tool. Use this for every mathematician and reviewer invocation throughout Phases 1-3.

• Fallback: Inline Reference Injection: references/agent-prompt-templates.md Section 2.

Apply only when an agent responds "TOOL_UNAVAILABLE: Read" -- re-invoke with reference content embedded directly in the prompt.

• Perspective + Synthesis Agent Templates (FULL mode swarm): Sections 3-4

of the same reference file. The decomposed swarm pipeline in Section 8b drives these.

---

7. Pre-Flight Validation Protocol -- Phase 0

MANDATORY: Before executing ANY workflow step, run pre-flight validation.

Step 0.1: Environment Validation

1. Read references/environment-setup.md Section 4 (Pre-Flight Validation Script) 2. Generate preflight_check.py and execute it via Bash 3. Interpret results per decision table:

Critical Rule: If pre-flight fails on dolfinx or petsc4py, do NOT attempt to generate simulation scripts. Help the user set up the environment first.

Step 0.2: Agent Availability Checks

Verify specialist agents are deployed:

1. Check cfd-mathematician SKILL.md exists at deployed path 2. Check cfd-reviewer SKILL.md exists at deployed path 3. [Removed] brainstorming-pm is no longer invoked as a sub-agent. FULL mode perspective agents are spawned directly by the orchestrator. 4. If agents missing: "New agent skills not synced. Run: ./sync-config.py push" 5. If agents remain unavailable: degrade gracefully (see Section 1)

Step 0.3: Reference File Checks

Verify all reference files exist:

• references/environment-setup.md
• references/physics-models.md
• references/mesh-generation-guide.md
• references/fenicsx-patterns.md
• references/validation-benchmarks.md
• references/troubleshooting-guide.md
• references/orchestrator-handoff-schema.md
• references/agent-loading-guide.md
• references/swarm-framing-templates.md
• references/agent-prompt-templates.md

Step 0.3a: Skill Base Path Resolution

Determine the absolute path to the cfd-bioreactor skill directory:

1. This SKILL.md is loaded from a known filesystem location. The directory containing this file is the skill base path. 2. Define: SKILL_BASE = "/Users/davidangelesalbores/repos/claude/claude-config/skills/cfd-bioreactor" 3. Store skill_base_path in the session state file (see Section 20). 4. All subsequent agent invocations include {skill_base_path}/references/ as the prefix for reference file paths. 5. All subsequent Read operations by the orchestrator use absolute paths: {skill_base_path}/references/{filename}

Step 0.4: Session Initialization

1. Create session directory: /tmp/cfd-bioreactor-session-{YYYYMMDD-HHMMSS}-{PID}/ 2. Create subdirectories: handoffs/, scripts/, logs/, perspectives/ 3. Initialize state file (see Section 20 for schema)

Step 0.5: Problem Parsing and Mode Selection

1. Parse user request: extract geometry type, fluid properties, species parameters, boundary conditions, target tier 2. Compute initial dimensionless estimates (Re, Pe, Da) if parameters available 3. Select mode based on tier + complexity heuristics (Section 5) 4. Determine run_purpose based on tier and session history:

• Tier 1-2: Always debug
• Tier 3-4: Check for prior completed sessions with same geometry hash.

If no prior session: debug. If prior debug session completed successfully: production (user can override).

• If resuming a session from a state file that lacks run_purpose, default

to debug.

• Store run_purpose in state file (see Section 20)

5. Apply run_purpose/mode precedence (auto-adjust conflicts):

Display any auto-adjustments to the user with rationale. 6. Present mode display to user (Section 5e), including run purpose:

   Problem: [description] (Tier [N])
   Available modes: DIRECT (10 min) | LITE (20 min) | FULL (40 min)
   Recommended: [MODE] ([reason])
   Run purpose: [DEBUG | PRODUCTION] ([reason])
   Override: You can select a different mode or run purpose.

Quality Gate 0: Pre-flight passes (dolfinx and petsc4py available). Session directory created. State file initialized. Mode selected. Run purpose determined.

---

8. Phase 1 -- Mesh Planning

Step 1.1: Swarm Discussion (FULL Mode Only)

1. Read references/swarm-framing-templates.md Swarm 1 template 2. Fill in placeholders with problem parameters (geometry, physical groups, memory budget) 3. Execute the decomposed swarm pipeline (Section 8b): a. Spawn 3 perspective agents (Mesh Engineer, Numerical Analyst, Computational Pragmatist) with the filled Swarm 1 challenge template b. Collect perspective outputs from {session_dir}/perspectives/phase1-*.md c. Spawn synthesis agent to combine perspectives into swarm-synthesis handoff 4. Check swarm output against quality threshold:

• At least 2 specific numerical recommendations
• At least 1 concrete alternative approach

5. If below threshold: log "Swarm 1 output below quality threshold; proceeding with mathematician analysis only"

Step 1.2: Mathematical Analysis

Invoke cfd-mathematician via Task tool:

• Input: Geometry description, target error tolerance, element type options,

swarm synthesis (if FULL mode)

• Loading instructions: Include per agent-loading-guide.md (physics-models.md

Sections 1-3, 6; validation-benchmarks.md Benchmark 1)

• Task: Recommend element type/order, estimate required mesh density, analyze

boundary layer resolution needs

• Timeout: 8 minutes
• On timeout: Use reference file defaults (linear triangles/tetrahedra, h estimated

from geometry scale)

Step 1.3: Adversarial Review

Invoke cfd-reviewer via Task tool:

• Input: Mesh plan from Steps 1.1-1.2, mathematician handoff output
• Loading instructions: Include per agent-loading-guide.md
• Task: Challenge mesh quality thresholds, refinement adequacy, memory feasibility,

STEP unit assumptions

• Timeout: 8 minutes
• On timeout: Proceed with APPROVED_WITH_WARNINGS, log "NO ADVERSARIAL REVIEW"

Quality Gate 1

Mechanical checks:

• Handoff YAML exists with required fields
• Memory estimate < 70% available RAM
• Element count > 0
• Mesh quality threshold (min scaled Jacobian > 0.1) is specified

Agent assessment:

• Reviewer approval_status is APPROVED or APPROVED_WITH_WARNINGS
• If APPROVED_WITH_WARNINGS: no CRITICAL-severity challenges
• If REJECTED: pass reviewer blocking_issues as HARD CONSTRAINTS to mathematician

(retry Step 1.2)

• Max 1 retry of mathematician after rejection. If still REJECTED: escalate to user

(see Section 16)

Handoff: Write Phase 1 mesh-plan handoff to {session_dir}/handoffs/phase1-mesh-plan.yaml

---

8b. Decomposed Swarm Pipeline (FULL Mode Only)

This section describes the two-phase decomposed pipeline that replaces the brainstorming-pm invocation for FULL mode swarm discussions. The orchestrator spawns perspective agents directly (respecting the Claude Code depth-1 constraint) and then spawns a synthesis agent to combine their outputs.

Why not brainstorming-pm? The brainstorming-pm skill is designed as a standalone orchestrator that spawns 5 parallel perspective agents via Task tool. When invoked as a sub-agent by cfd-bioreactor, brainstorming-pm CANNOT use the Task tool (Claude Code platform constraint: sub-agents cannot spawn sub-agents). This would reduce its multi-perspective pipeline to single-context sequential operation, losing the diversity that makes swarm discussion valuable.

Architecture:

cfd-bioreactor (depth 0)
  |-- Task --> Perspective Agent: Numerical Analyst   (depth 1)
  |-- Task --> Perspective Agent: Mesh Engineer        (depth 1)
  |-- Task --> Perspective Agent: Physical Modeler     (depth 1)
  |-- Task --> Perspective Agent: Computational Pragmatist (depth 1)
  |-- Task --> Perspective Agent: Validation Strategist (depth 1, optional)
  |
  |-- [collect perspective outputs from files]
  |
  |-- Task --> Synthesis Agent                         (depth 1)

Phase D1: Spawn Perspective Agents

1. Read references/swarm-framing-templates.md for the appropriate swarm template (Swarm 1, 2, or 3) and fill in placeholders with problem parameters 2. Spawn 3-5 perspective agents as parallel Task calls. Fill the verbatim Perspective Agent Task Template in references/agent-prompt-templates.md Section 3 (one filled copy per agent, with its {perspective_name} and the challenge text from references/swarm-framing-templates.md) and pass each to the Task tool. 3. Apply per-agent timeout: 5 minutes 4. Minimum agents required: 3 of 5 5. If fewer than 3 complete: log warning, skip swarm synthesis, proceed as LITE mode

Perspective definitions (see swarm-framing-templates.md Section 6 for full prompts):

For Phase 1 (Mesh): prioritize Mesh Engineer, Numerical Analyst, Computational Pragmatist (3 agents minimum) For Phase 2 (Flow): prioritize Numerical Analyst, Physical Modeler, Computational Pragmatist For Phase 3 (Transport): prioritize Numerical Analyst, Physical Modeler, Validation Strategist

Phase D2: Synthesis Agent

After all perspective agents complete (or timeout):

1. Collect perspective outputs from {session_dir}/perspectives/phase{N}-*.md 2. Spawn a synthesis agent via Task tool. Fill the verbatim Synthesis Agent Task Template in references/agent-prompt-templates.md Section 4 (it includes the exact swarm-synthesis handoff YAML format the agent must emit) and pass it to the Task tool. 3. Apply synthesis agent timeout: 5 minutes 4. If synthesis agent fails or times out: extract insights manually from perspective files (orchestrator reads files and summarizes)

Swarm Quality Check

After Phase D2 completes, validate the swarm synthesis:

• At least 2 specific numerical recommendations in convergent_insights
• At least 1 concrete alternative in divergent_alternatives
• If below threshold: log "Swarm output for Phase {N} below quality threshold"

and proceed without swarm input

---

9. Phase 2 -- Flow Solver Planning, Code Generation, and Execution

Step 2.1: Swarm Discussion (FULL Mode Only)

Same pattern as Step 1.1 with Swarm 2 template from references/swarm-framing-templates.md. Perspective agents for Phase 2: Numerical Analyst, Physical Modeler, Computational Pragmatist (see Section 8b phase-perspective mapping).

Step 2.2: Mathematical Analysis

Invoke cfd-mathematician via Task tool:

• Input: Flow parameters (Re, fluid properties), mesh plan from Phase 1,

swarm synthesis (if FULL mode)

• Task: Write variational formulation (weak form), verify Taylor-Hood inf-sup

stability, estimate convergence order, recommend continuation strategy if Re > 10

• Timeout: 8 minutes

Step 2.3: Adversarial Review

Invoke cfd-reviewer via Task tool:

• Input: Flow plan from Steps 2.1-2.2, mathematician handoff output
• Task: Challenge BC physical consistency, solver convergence criteria, check for

missing pressure reference point, verify Newton vs. Picard choice, check memory estimates

• Timeout: 8 minutes

Quality Gate 2a: Flow Plan Approved

Same pattern as Quality Gate 1:

• Mechanical checks + reviewer approval
• Max 1 retry on rejection, then escalate to user

Step 2.4: Code Generation (Orchestrator-Owned)

1. Read reference file sections per agent-loading-guide.md Phase 1 + Phase 2 maps 2. Generate mesh script from Phase 1 validated plan:

• Apply mesh sizing defaults from Section 15b based on run_purpose:

debug runs use the coarse element count range for the current tier; production runs use the fine element count range

• Geometric feature resolution check: ensure >= 3 elements across the thinnest

geometric feature (membranes, thin channels). If the debug element count range cannot achieve this, increase to the minimum that satisfies this constraint and warn the user.

• Import STEP geometry via gmsh.model.occ.importShapes() (with error handling)

OR construct parametric geometry via gmsh built-in/OCC kernel

• Call gmsh.model.occ.synchronize() after ALL OCC operations
• Define physical groups (inlet, outlet, walls, membrane, cell_region)
• Apply mesh refinement near walls and membrane interfaces
• Estimate memory requirements before meshing
• Check mesh quality after generation
• Convert to DOLFINx via dolfinx.io.gmsh.model_to_mesh()

3. Generate flow solver script from Phase 2 validated plan:

• Define Taylor-Hood P2/P1 function space
• Apply boundary conditions from physical groups
• For Re < 1: assemble Stokes variational form (linear solve)
• For Re >= 1: Navier-Stokes with Newton continuation:

a. Solve Stokes for initial guess b. Ramp Re through intermediate values if Re > 10 c. Use previous solution as initial guess for each step

• Select solver: MUMPS for < 50K DOFs, GMRES+ILU otherwise
• Include solver progress monitor for 3D
• Save velocity + pressure fields to XDMF checkpoint

4. Apply ALL code generation rules from Section 13 (Code Generation Protocol) 5. Run syntax check: python3 -c "import ast; ast.parse(open('script.py').read())"

Step 2.5: Execution and Validation

1. Execute mesh + flow scripts (Bash for 2D, Write for 3D) 2. Validate:

• Solver converged (residual < tolerance)
• Mass conservation: |integral u.n ds_inlet + integral u.n ds_outlet| < 1e-6
• No NaN or Inf in solution arrays
• No unphysical flow patterns (reversed flow at inlet, etc.)
• If Poiseuille comparison available: L2 error < 1% (P2 elements give ~1e-12 for

quadratic solutions; see references/validation-benchmarks.md Benchmark 1)

Step 2.5b: Self-Correction Loop (If Execution or Validation Fails)

1. Collect error output (traceback, log messages, solver convergence history) 2. Classify error type: import_error, mesh_error, solver_divergence, assertion_failure, numerical_instability, wall_clock_timeout 3. Simple errors (import path, syntax, assertion): Fix directly and retry (max 2 direct retries) 4. Complex errors (solver divergence, numerical instability): Invoke cfd-reviewer in error diagnosis mode:

• Pass: error output, error_history from previous attempts, mathematical specification
• Reviewer classifies root cause and recommends specific fix
• Check error_history to avoid recommending previously-failed fixes

4b. Wall-clock timeout: If the simulation hit the kill threshold from Section 15b:

• If run_purpose == debug: The mesh is likely too fine. Reduce element count

by 50% and retry. Do NOT escalate to cfd-reviewer for this error type.

• If run_purpose == production: Invoke cfd-reviewer for diagnosis. Common

causes: insufficient Newton continuation steps, solver parameter tuning needed. 5. Regenerate code with reviewer feedback applied 6. Track error_history (append each attempt: error_type, error_message, fix_attempted, fix_outcome) 7. After 3 total retries: Escalate to user with full error summary:

   EXECUTION FAILED after 3 attempts.
   Error type: [classification]
   Attempts:
   1. [fix attempted] -> [outcome]
   2. [fix attempted] -> [outcome]
   3. [fix attempted] -> [outcome]
   Recommendation: [suggested next step]

Quality Gate 2b: Flow validation passes (mass conservation, convergence, no NaN/Inf).

Handoff: Write Phase 2 flow-result handoff to {session_dir}/handoffs/phase2-flow-result.yaml

---

10. Phase 3 -- Transport Planning, Code Generation, and Execution

Mirrors Phase 2 structure with transport-specific content.

Step 3.1: Swarm Discussion (FULL Mode Only)

Same pattern as Step 1.1 with Swarm 3 template from references/swarm-framing-templates.md. Perspective agents for Phase 3: Numerical Analyst, Physical Modeler, Validation Strategist (see Section 8b phase-perspective mapping).

Step 3.2: Mathematical Analysis

Invoke cfd-mathematician via Task tool:

• Input: Transport parameters (Pe, Da, species properties), velocity field summary

from Phase 2, swarm synthesis (if FULL mode)

• Task: Write transport variational form with SUPG, analyze stabilization parameter

(Pe-dependent xi formula), verify regularized MM convergence properties, estimate transport convergence order

• Timeout: 8 minutes

Step 3.3: Adversarial Review

Invoke cfd-reviewer via Task tool:

• Input: Transport plan from Steps 3.1-3.2, mathematician handoff output
• Task: Challenge SUPG parameter for overflow risk (Pe > 710 with cosh/sinh), check

MM regularization epsilon adequacy, verify Robin BC formulation for membrane, check Newton solver convergence expectations

• Timeout: 8 minutes

Quality Gate 3a: Transport Plan Approved

Same pattern as Quality Gates 1 and 2a:

• Mechanical checks + reviewer approval
• Max 1 retry on rejection, then escalate to user

Step 3.4: Code Generation (Orchestrator-Owned)

1. Read reference file sections per agent-loading-guide.md Phase 3 map 2. Apply mesh sizing from Section 15b for any transport-specific mesh refinement (e.g., near-membrane refinement for SUPG accuracy). Use run_purpose to select the appropriate element count range. Ensure >= 3 elements across thin features (membranes, boundary layers) even in debug mode. 3. Generate transport script from Phase 3 validated plan:

• Define P1 function space for concentration
• Estimate Peclet number: Pe = |u| * h / (2D)
• Implement SUPG stabilization (see Section 13 -- ALWAYS rules)
• Implement regularized Michaelis-Menten (see Section 13 -- ALWAYS rules)
• Implement membrane permeation as Robin BC (Fick's law)
• Use Newton solver (nonlinear due to Michaelis-Menten)
• Import NonlinearProblem from dolfinx.fem.petsc, NewtonSolver from dolfinx.nls.petsc
• Include post-solve checks: negative concentration warning, species conservation

4. Apply ALL code generation rules from Section 13 5. Run syntax check

Step 3.5: Execution and Validation

1. Execute transport script (Bash for 2D, Write for 3D) 2. Validate:

• Newton solver converged
• Species conservation: |inlet_flux + outlet_flux + membrane_flux + total_reaction| < 1%

of total reaction

• min(c) check: warn if significantly negative (suggests insufficient mesh resolution)

Step 3.5b: Self-Correction Loop

Same pattern as Step 2.5b, including the wall_clock_timeout error type. Invoke cfd-reviewer in error diagnosis mode for complex errors. For wall-clock timeouts: if run_purpose == debug, reduce element count by 50%; if run_purpose == production, invoke cfd-reviewer. Max 3 total retries, then escalate to user.

Quality Gate 3b: Transport validation passes (species conservation, Newton convergence, no large negative concentration regions).

Handoff: Write Phase 3 transport-result handoff to {session_dir}/handoffs/phase3-transport-result.yaml

---

11. Phase 4 -- Post-Processing and Visualization

Orchestrator-owned (no agent involvement).

Step 4.1: Visualization Scripts

1. Read references/fenicsx-patterns.md Section 12 (PyVista patterns) 2. Generate visualization script with:

• Velocity magnitude contour plot
• O2 concentration field with depletion zones highlighted
• Streamlines from inlet (3D)
• Cross-section slices at multiple positions (3D)
• Centerline line probe (O2 vs. position)
• Screenshots saved to PNG files
• VTK files saved for ParaView
• Headless fallback: if not pyvista.system_supports_plotting(): pyvista.OFF_SCREEN = True

3. Execute 2D visualizations; save 3D scripts for manual execution

Step 4.2: Mesh Convergence Study (Optional -- for Publication Quality)

1. Read references/validation-benchmarks.md Section 4 (convergence protocol) 2. Generate convergence study script:

• 4 mesh refinement levels: h, h/2, h/4, h/8 (or 3 for 3D)
• Solve on each level
• Compute QoI: outlet average O2, max velocity, min O2 in cell region
• Richardson extrapolation for grid-independent value
• Convergence rate: p = log(e_{i-1}/e_i) / log(2)
• Non-monotone convergence detection and warning

Quality Gate 4a (convergence study):

• Monotone convergence (error decreases with refinement)
• Convergence rate matches expected order (O(h^2) for P1 concentration, O(h^3) for P2 velocity)
• QoI changes < 1% between last two levels

---

12. Phase 5 -- Summary and Deliverables

1. Collect all handoff documents from {session_dir}/handoffs/ 2. Report to user:

• What was simulated (geometry, physics, parameters)
• Key agent decisions (mathematician recommendations, reviewer concerns, resolutions)
• Validation results (mass conservation, species conservation, convergence)
• Files produced (scripts, XDMF checkpoints, VTK files, PNG images)

3. Offer to preserve or clean up session directory 4. Update state file with status: "complete"

---

13. Code Generation Protocol

When generating FEniCSx Python scripts, ALWAYS follow these rules. These rules are embedded inline because they guard against code-level errors that produce silently wrong results. They must always be in the orchestrator's context.

ALWAYS Include

1. Version assertion header at the top of every script:

   import importlib.metadata as _meta
   _ver = tuple(int(x) for x in _meta.version("fenics-dolfinx").split(".")[:2])
   assert _ver >= (0, 10), f"Requires FEniCSx >= 0.10. See references/environment-setup.md."

2. Reproducibility metadata header: date, dolfinx version, basix version, gmsh version, numpy version, Python version, OS, all physical parameters

3. Complete, self-contained scripts: Every script must run standalone from the command line. Never generate code fragments that require manual assembly.

4. Inline comments explaining each step: what it does, why it is needed, what the expected output is

5. Regularized Michaelis-Menten: Always use c_pos = (c + sqrt(c^2 + eps^2)) / 2 with eps = 1e-10 * c_inlet. Never use the raw c / (Km + c) form.

6. SUPG stabilization for all species transport (check Pe first, but enable by default): Use the numerically stable form xi = conditional(gt(Pe, 1.0), 1.0 - 1.0/Pe, Pe/3.0). Never use cosh/sinh (overflows for Pe > 710).

7. Post-solve quality checks: mass conservation, species conservation, negative concentration warning, NaN/Inf detection

8. Save results to XDMF (for checkpointing) and VTK/PVD (for PyVista/ParaView):

   with io.XDMFFile(comm, "results.xdmf", "w") as xdmf:
       xdmf.write_mesh(domain)
       xdmf.write_function(uh, 0.0)

NEVER Do

• Generate code fragments (always complete scripts)
• Use from dolfin import * (legacy FEniCS, not FEniCSx)
• Use dolfinx.io.gmshio (v0.10 renamed to dolfinx.io.gmsh)
• Use raw Michaelis-Menten c / (Km + c) without regularization
• Use cosh(Pe)/sinh(Pe) for SUPG parameter (numerical overflow)
• Skip physical group definition before meshing
• Forget gmsh.model.occ.synchronize() after OCC operations
• Use equal-order elements (P1/P1) for flow without stabilization
• Omit the pressure reference point (at least one pressure BC needed)

Pattern Library

Use code patterns from references/fenicsx-patterns.md as the primary source for all generated code. These patterns have been validated against FEniCSx v0.10 API.

Use parameter values from references/physics-models.md lookup tables for bioreactor-specific constants (fluid properties, diffusion coefficients, Michaelis-Menten kinetics by cell type).

Additional Code Generation Rules

• Taylor-Hood P2/P1 for Stokes/Navier-Stokes flow (standard inf-sup stable pair)
• Newton continuation for Re > 10: Solve Stokes for initial guess, ramp Re through

intermediate values [1, 10, 50, target], use previous solution as initial guess

• MUMPS for < 50K DOFs (direct solver), GMRES+ILU for larger problems (iterative)
• NonlinearProblem from dolfinx.fem.petsc, NewtonSolver from dolfinx.nls.petsc
• *eps = 1e-10 c_inlet** for Michaelis-Menten regularization parameter

---

14. Error Handling Protocol

Simulation Error Handling

When errors occur during simulation, follow this diagnostic sequence:

Solver Divergence Recovery Sequence

If the Newton solver diverges, try these fixes in order:

1. Use Stokes solution as initial guess (most common fix) 2. Ramp Reynolds number through intermediate values: [1, 10, 50, target] 3. Reduce relaxation factor to 0.5 (then 0.3, then 0.1) 4. Switch to Picard iteration instead of Newton 5. Refine mesh near boundaries and high-gradient regions 6. If all fail: suggest the user simplify the geometry or reduce Re

Always provide actionable diagnostic messages. Never just report an error code.

Orchestration Error Handling

---

15. Timeout Configuration

Per-Agent Timeouts

Per-Phase Execution Timeouts

Per-Phase Agent Invocation Budgets (FULL Mode)

These budgets include all agent invocations within the phase (perspectives, mathematician, reviewer, and any retries). If the budget is exceeded before reviewer approval, proceed with the current best plan and log: "Phase {N} time budget exceeded. Proceeding with current plan."

Cumulative budget: If total agent invocation time across Phases 1-3 exceeds 90 minutes in FULL mode, auto-downgrade remaining phases to LITE mode (skip perspective agents) and summarize all previous handoffs to 200-token summaries.

Rule: If estimated runtime exceeds 10 minutes, generate the script and instruct the user to run it manually. Do not attempt to execute long-running simulations via Bash.

For 3D simulations: Always include MPI instructions:

mpirun -np 4 python simulation.py

---

15b. Simulation Runtime Management

15b.1 Debug vs. Production Mesh Sizing

Simulation mesh sizing defaults are tied to the complexity tier and run_purpose. The orchestrator determines run_purpose during Phase 0 (Step 0.5) based on user intent and tier selection. There are two run purpose levels:

• debug: First run, error recovery, parameter exploration. Use coarsest mesh that

captures qualitative physics.

• production: Publication-quality results. Use fine mesh with convergence study.

Tier-aware mesh sizing defaults:

Geometric feature resolution rule: Regardless of run_purpose, ensure >= 3 elements across the thinnest geometric feature (membranes, thin channels, boundary layers). If the default element count range for a given tier and run_purpose cannot achieve this, increase the element count to the minimum that satisfies this constraint and warn the user: "Debug mesh increased from [default] to [adjusted] elements to resolve [feature] (>= 3 elements across thinnest feature required)."

Orchestrator behavior: 1. During Phase 0 (Step 0.5), set run_purpose based on tier and session history 2. During Phase 1 (mesh planning), pass run_purpose to cfd-mathematician with the corresponding element count range as a constraint 3. During code generation (Steps 2.4/3.4), use the mesh sizing defaults from this table unless the mathematician or reviewer recommends otherwise 4. If run_purpose == debug, skip Phase 4 mesh convergence study entirely

Override: The user can override run purpose at any time:

Current run purpose: debug (Tier 3, first run)
Override to production? This will increase mesh from ~10K to ~100K elements
and runtime from ~2 minutes to ~30-60 minutes.

15b.2 Solver Convergence Assessment

During execution (Steps 2.5, 3.5), monitor solver convergence and apply these diagnostic thresholds:

Newton solver convergence:

Krylov solver convergence (GMRES/CG inner solves):

Action on convergence signals:

• Healthy: Continue normally
• Warning: Log warning, continue, but flag for user in Phase 5 summary.

If run_purpose == debug, suggest "convergence marginal on debug mesh -- may improve on finer production mesh"

• Diverging: Enter self-correction loop (Step 2.5b/3.5b). If first attempt

on a debug mesh, do NOT refine mesh as first fix -- instead try solver parameter adjustments (relaxation, Picard, continuation ramp)

15b.3 Mesh Convergence Decision Guide

Phase 4 mesh convergence study (Section 11, Step 4.2) is currently marked as optional. This section provides guidance on WHEN to run it.

Run mesh convergence study when:

• run_purpose == production AND tier >= 3
• User explicitly requests publication-quality results
• QoI sensitivity: if a 10% mesh coarsening changes QoI by > 5%, convergence

study is needed

Skip mesh convergence study when:

• run_purpose == debug (always skip)
• Tier 1-2 (analytical solution available for validation instead)
• User explicitly declines ("I just need qualitative results")

Convergence assessment criteria (extending Quality Gate 4a):

• QoI change between last two refinement levels < 1% (existing criterion)
• Observed convergence rate within 0.5 of theoretical order (new criterion):
• P1 concentration: expect p ~ 2.0, accept 1.5-2.5
• P2 velocity: expect p ~ 3.0, accept 2.5-3.5
• If convergence rate is anomalously low (p < 1.0): flag as potential

singularity or insufficient boundary layer resolution

15b.4 Wall-Clock Time Limits

Prevent simulations from running indefinitely by setting wall-clock expectations per tier and detecting stalls.

Expected wall-clock times (serial execution, single core; times are for simulation execution only, not including agent discussion time):

Kill threshold: If wall-clock time exceeds 3x the expected total for the tier and run purpose, the simulation is likely stalled or diverging.

Orchestrator behavior when kill threshold is reached: 1. If running via Bash: timeout will have already killed the process 2. If running manually (3D production): include a timeout check in the generated script:

   import time
   _start_time = time.time()
   _kill_threshold_s = {kill_threshold_seconds}

   # Inside solver loop:
   if time.time() - _start_time > _kill_threshold_s:
       print(f"TIMEOUT: Simulation exceeded {_kill_threshold_s}s kill threshold.")
       print("Likely causes: mesh too fine for debug, solver diverging, or")
       print("insufficient Newton continuation steps.")
       print("Try: reduce mesh size, increase continuation ramp steps,")
       print("or switch to Picard iteration.")
       sys.exit(1)

3. Report timeout cause in Phase 5 summary with specific recommendations

Stall detection heuristic (for solver progress monitoring):

• If solver has not reduced residual by > 10x in the last 5 minutes of wall-clock

time (for production runs) or 30 seconds (for debug runs), flag as stalled

• Stall is different from slow convergence: stall means no progress at all

---

16. Quality Gate Escalation Protocol

When the mathematician and reviewer cannot agree after 1 retry, present the disagreement to the user in this structured format:

AGENT DISAGREEMENT: [Phase Name]

Mathematician recommends: [1-2 sentence summary]
Rationale: [1 sentence]

Reviewer objects: [1-2 sentence summary]
Concern: [1 sentence]

Options:
(A) Proceed with mathematician's recommendation (risk: [reviewer's concern])
(B) Proceed with reviewer's preferred approach (trade-off: [mathematician's concern])
(C) Use conservative approach: [orchestrator's auto-generated compromise]
(D) Provide your own parameters

Recommended: (C) [brief justification]

Conservative compromise rule: Always take the more cautious choice from both agents -- finer mesh, more stable solver, lower tolerance, more regularization.

Escalation sequence: 1. First rejection: Reviewer's blocking_issues become HARD CONSTRAINTS for mathematician retry 2. After 1 retry still REJECTED: Escalate to user with the structured format above 3. User selects option (A/B/C/D) or provides custom parameters 4. Orchestrator proceeds with user's choice

---

17. Orchestrator Context Management

The orchestrator loads reference file sections on a per-phase basis to control context window consumption. Agent discussion details stay in handoff files on disk -- they are NOT loaded back into the orchestrator context.

Previous phase summaries: Max 200 tokens per phase. Include only: key decisions, approval status, validation metrics. Not full discussion content.

Section 13 NEVER rules are always in context (they are embedded inline in this SKILL.md, not in a reference file).

---

18. Phase Dependencies

Rule: If Phase N is re-executed, all phases > N are invalidated and must be re-run.

---

19. Complexity Tiers

Tier 1: Validation (5-15 minutes) --> DIRECT Mode

2D Poiseuille flow. Tests the entire pipeline without any user-specific physics. Always start here when the user has never run a FEniCSx simulation before. Default mesh sizing: 200-500 elements (2D). Always debug purpose.

See: examples/01-2d-channel-flow.md

Tier 2: Single Physics (15-30 minutes) --> LITE Mode

2D user geometry + O2 transport with Michaelis-Menten and membrane BCs. Demonstrates the full two-phase workflow on a simple geometry. Default mesh sizing: 500-2,000 elements (2D). Always debug purpose.

See: examples/02-2d-oxygen-transport.md

Tier 3: Coupled Multiphysics (1-4 hours) --> FULL Mode

2D or 3D STEP geometry + full coupled physics. Production-quality simulation with iterative solvers, Newton continuation, and validation suite. Default mesh sizing: Debug: 5,000-20,000 elements (3D). Production: 50,000-200,000 elements (3D). First run always defaults to debug. Override to production after validation passes.

See: examples/03-3d-cartridge-template.md

Tier 4: Production 3D (4-24 hours) --> FULL Mode

Full 3D cartridge with mesh convergence study, publication-quality visualization, and complete reproducibility metadata. Always run with MPI on a workstation. Default mesh sizing: Debug: 10,000-50,000 elements (3D). Production: 200,000-1,000,000+ elements (3D). First run always defaults to debug. Override to production after validation passes.

Progression rule: Suggest the user start at Tier 1 and progress upward. Do not jump to Tier 3 or 4 without first validating the environment at Tier 1.

Mode mapping:

---

20. Session State File Schema

State file location: /tmp/cfd-bioreactor-state-{session-id}.yaml

session:
  session_id: "20260221-143000-12345"
  session_dir: "/tmp/cfd-bioreactor-session-20260221-143000-12345"
  skill_base_path: "/Users/davidangelesalbores/repos/claude/claude-config/skills/cfd-bioreactor"
  mode: "LITE"                    # DIRECT | LITE | FULL
  tier: 2                         # 1-4
  run_purpose: "debug"            # debug | production
  handoff_declined: false         # Set to true if user declines programming-pm handoff
  current_phase: 2                # 0-5
  status: "running"               # running | paused | complete | failed

phases_completed: [0, 1]          # List of completed phase numbers
phases_valid: [0, 1]              # List of phases whose results are still valid

phase_results:
  phase_0:
    status: "complete"
    handoff_path: null
    timestamp: "2026-02-21T14:30:00Z"
  phase_1:
    status: "complete"
    handoff_path: "{session_dir}/handoffs/phase1-mesh-plan.yaml"
    reviewer_status: "APPROVED_WITH_WARNINGS"
    retries_used: 0
    timestamp: "2026-02-21T14:35:00Z"
  phase_2:
    status: "running"
    handoff_path: null
    reviewer_status: null
    retries_used: 0
    execution_retries_used: 0     # For self-correction loop

timestamps:
  started: "2026-02-21T14:30:00Z"
  last_updated: "2026-02-21T14:40:00Z"

errors: []                        # List of non-fatal errors/warnings encountered

---

21. Session Resume Protocol

On startup, check for existing state files:

1. Search for /tmp/cfd-bioreactor-state-*.yaml with status != "complete" 2. If found: display last completed phase, offer to resume

   Found incomplete session from [timestamp].
   Last completed phase: Phase [N] ([phase_name])
   Mode: [DIRECT/LITE/FULL], Tier: [N]
   Resume from Phase [N+1]? (yes/no/restart)

3. If user chooses resume:

• Phase 0: Always re-run (fast, catches environment changes)
• Phases 1-4: Resumable if handoff files exist at recorded paths
• Verify handoff file existence before resuming

4. If user chooses restart: Create new session, archive old state file

---

22. Integration with Other Skills

22b. Programming-PM Handoff Criteria

The cfd-bioreactor orchestrator generates standalone simulation scripts. When the user's needs evolve beyond one-off scripts into maintained software, recommend handoff to programming-pm for software development coordination.

Handoff Trigger Conditions

Recommend programming-pm when ANY of the following are detected:

Context Filtering Rules

To prevent false-positive handoff triggers, apply these filters:

1. Phase filtering: NEVER evaluate handoff triggers during Phases 2-4 (active computation). Only evaluate at phase boundaries (Phase 0, Phase 5) or during user interaction between phases.

2. Compound requirement: Require keywords from 2+ different trigger categories (e.g., "library" + "test") OR a single unambiguous keyword that has no alternate interpretation in CFD context (e.g., "package this as a library", "set up CI").

3. Exclusion patterns: Do NOT treat these as handoff triggers:

• "test" in "test this simulation", "test run", "test different parameters"

(means "try", not "unit test")

• "git" in "let me check" or "let me get" (not version control)
• "version" in "FEniCSx version", "Python version" (not software versioning)
• "module" in "Python module not found" (import error, not packaging)

4. Handoff decline suppression: If user selects option (B) Continue in the handoff recommendation, set handoff_declined: true in session state. Do NOT re-suggest programming-pm handoff for the remainder of this session.

Handoff Recommendation Format

When a trigger condition passes context filtering, present to the user:

HANDOFF RECOMMENDATION: Software Development Coordination

Your simulation work is evolving into a software project:
- Trigger: [which trigger was detected]
- Current state: [N scripts generated, current tier/mode]

Recommendation: Hand off to programming-pm for:
- Architecture design (module structure, shared utilities)
- Testing strategy (unit tests for simulation components)
- Code review and quality gates
- Version control integration

Options:
(A) Hand off to programming-pm now (recommended)
(B) Continue with cfd-bioreactor (standalone scripts only)
(C) Complete current simulation phase, then hand off

What Gets Handed Off

When handing off to programming-pm, provide:

• All generated scripts from {session_dir}/scripts/
• Session state summary (tier, mode, run_purpose, phases completed)
• Physical parameters and validation results
• List of scripts and their purposes

When NOT to Hand Off

Do NOT recommend programming-pm for:

• Standard tier progression (Tier 1 -> 2 -> 3 -> 4)
• Mesh convergence studies (these are cfd-bioreactor Phase 4)
• One-off parameter variations (just re-run with different inputs)
• Visualization enhancements (cfd-bioreactor Phase 4)
• Any session where handoff_declined is already true

---

23. Reference Files

This skill includes the following reference files. Read them as needed during workflow execution -- do not load all at once. See references/agent-loading-guide.md for section-level loading maps per agent.

---

24. Notes

• Version pinning: All code patterns target FEniCSx v0.10. If a newer version is

released, the version assertion will catch the mismatch. Update patterns before using with a new API version.

• SI units everywhere: All parameters are in SI (kg, m, s, mol, Pa) unless the user

explicitly specifies otherwise. If a STEP file appears to be in mm (bounding box > 1 m), warn the user about potential unit mismatch.

• Serial execution by default: Generated scripts run in serial (single core). For 3D

production runs, add MPI instructions. The scripts are MPI-aware (use MPI.COMM_WORLD throughout) and work correctly with mpirun -np N.

• Reproducibility: Every generated script includes a metadata header with all version

numbers, parameters, and settings. Results can be reproduced by re-running the same script in the same environment.

• This skill generates code but does not modify system packages: It writes Python

scripts to the user's working directory. It never installs packages, modifies micromamba environments, or changes system configuration.

• Physical group convention: inlet=1, outlet=2, walls=3, membrane=4, cell_region=5,

fluid_volume=10. Consistent across all examples and generated scripts.

• gmsh synchronize: After ANY gmsh.model.occ.* operation (addBox, addCylinder,

importShapes, etc.), ALWAYS call gmsh.model.occ.synchronize() before accessing entities, assigning physical groups, or meshing. Forgetting this is the single most common gmsh error.

• v2.0 multi-agent orchestrator: This is v2.0 of the cfd-bioreactor skill. v1.0 was

a single-context generate-explain-validate skill; v2.0 is a multi-agent orchestrator with cfd-mathematician, cfd-reviewer, and decomposed perspective agent coordination.

Example 01: 2D Poiseuille Channel Flow (Tier 1 Validation)

Overview

What: Fully developed laminar flow (Poiseuille flow) between two parallel plates in a 2D rectangular channel. This is the simplest possible CFD validation case.

Why: The Poiseuille solution is a quadratic polynomial in the transverse coordinate. P2 (quadratic) velocity elements reproduce this exactly, so the velocity L2 error should be at machine precision (~1e-12). This validates that the entire pipeline -- geometry, mesh, function spaces, boundary conditions, solver, and post-processing -- is working correctly. Pressure convergence (P1 elements, O(h^2)) provides a meaningful convergence rate demonstration.

Time: 1-3 minutes on any workstation.

Complexity: Tier 1 (flow only, no species transport).

---

Complete Python Script

#!/usr/bin/env python3
"""
Poiseuille flow validation: 2D channel between parallel plates.

Analytical solution: u(y) = (dP / (2 mu L)) * y * (H - y)
P2 elements reproduce quadratic solutions exactly, so velocity error ~ 1e-12.
Pressure convergence (P1, O(h^2)) demonstrates mesh convergence methodology.

CFD Bioreactor Skill -- Tier 1 Validation
"""

# ── Version assertion ────────────────────────────────────────────────
import importlib.metadata as _meta

_dolfinx_ver = tuple(int(x) for x in _meta.version("fenics-dolfinx").split(".")[:2])
assert _dolfinx_ver >= (0, 10), (
    f"This script requires FEniCSx >= 0.10 (found {_meta.version('fenics-dolfinx')}). "
    "See references/environment-setup.md for installation instructions."
)

# ── Reproducibility header ───────────────────────────────────────────
import datetime, platform, dolfinx, basix, gmsh, pyvista, numpy as np

print("=" * 60)
print("Poiseuille Flow Validation")
print(f"Date:     {datetime.datetime.now().isoformat()}")
print(f"dolfinx:  {dolfinx.__version__}")
print(f"basix:    {basix.__version__}")
print(f"gmsh API: {gmsh.GMSH_API_VERSION}")
print(f"numpy:    {np.__version__}")
print(f"Python:   {platform.python_version()}")
print(f"OS:       {platform.platform()}")
print("=" * 60)

# ── Imports ──────────────────────────────────────────────────────────
from dolfinx import fem, mesh, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import gmsh as gmsh_io
from basix.ufl import element, mixed_element
import ufl
from mpi4py import MPI

# ── Physical parameters (SI units) ──────────────────────────────────
H = 1.0e-3       # channel height [m]
L = 10.0e-3      # channel length [m]
mu = 1.0e-3      # dynamic viscosity [Pa.s]
dP = 100.0       # pressure drop [Pa]

u_max = dP * H**2 / (8.0 * mu * L)
Re = 1000.0 * u_max * H / mu  # rho = 1000 kg/m^3
print(f"\nParameters: H={H*1e3:.1f} mm, L={L*1e3:.1f} mm, mu={mu:.1e} Pa.s, dP={dP:.0f} Pa")
print(f"u_max = {u_max:.4f} m/s, Re = {Re:.3f} (Stokes regime)")

# ── Geometry: parametric rectangle via gmsh built-in kernel ──────────
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)  # suppress terminal output
gmsh.model.add("channel")

# Rectangle: origin (0, 0), width L, height H
rect = gmsh.model.occ.addRectangle(0.0, 0.0, 0.0, L, H)
gmsh.model.occ.synchronize()

# Identify boundary lines by bounding box
eps = 1e-8 * min(H, L)
lines = gmsh.model.getEntities(dim=1)

inlet_lines, outlet_lines, bottom_lines, top_lines = [], [], [], []
for dim, tag in lines:
    xmin, ymin, _, xmax, ymax, _ = gmsh.model.getBoundingBox(dim, tag)
    if abs(xmin) < eps and abs(xmax) < eps:
        inlet_lines.append(tag)
    elif abs(xmin - L) < eps and abs(xmax - L) < eps:
        outlet_lines.append(tag)
    elif abs(ymin) < eps and abs(ymax) < eps:
        bottom_lines.append(tag)
    elif abs(ymin - H) < eps and abs(ymax - H) < eps:
        top_lines.append(tag)

# Physical groups (MANDATORY before meshing)
INLET = 1; OUTLET = 2; BOTTOM_WALL = 3; TOP_WALL = 4; FLUID = 10
gmsh.model.addPhysicalGroup(1, inlet_lines,  INLET,       name="inlet")
gmsh.model.addPhysicalGroup(1, outlet_lines,  OUTLET,      name="outlet")
gmsh.model.addPhysicalGroup(1, bottom_lines, BOTTOM_WALL, name="bottom_wall")
gmsh.model.addPhysicalGroup(1, top_lines,    TOP_WALL,    name="top_wall")
gmsh.model.addPhysicalGroup(2, [rect],       FLUID,       name="fluid")

# Mesh
N_across = 16  # elements across channel height
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", H / N_across)
gmsh.model.mesh.generate(2)
print(f"\nMesh: {len(gmsh.model.mesh.getNodes()[0]) // 3} nodes")

# ── Convert to DOLFINx ──────────────────────────────────────────────
mesh_data = gmsh_io.model_to_mesh(gmsh.model, MPI.COMM_WORLD, rank=0, gdim=2)
domain = mesh_data.mesh
cell_tags = mesh_data.cell_tags
facet_tags = mesh_data.facet_tags
gmsh.finalize()

# ── Taylor-Hood P2/P1 function space ────────────────────────────────
cell_name = domain.topology.cell_name()
P2 = element("Lagrange", cell_name, 2, shape=(2,))
P1 = element("Lagrange", cell_name, 1)
TH = mixed_element([P2, P1])
W = fem.functionspace(domain, TH)

# ── Boundary conditions ─────────────────────────────────────────────
fdim = domain.topology.dim - 1

# Analytical inlet velocity profile: u(y) = (dP/(2*mu*L)) * y * (H - y)
def inlet_profile(x):
    vals = np.zeros((2, x.shape[1]))
    vals[0] = (dP / (2.0 * mu * L)) * x[1] * (H - x[1])
    return vals

# Inlet: prescribed parabolic velocity
V_sub, _ = W.sub(0).collapse()
u_inlet = fem.Function(V_sub)
u_inlet.interpolate(inlet_profile)
inlet_facets = facet_tags.find(INLET)
inlet_dofs = fem.locate_dofs_topological((W.sub(0), V_sub), fdim, inlet_facets)
bc_inlet = fem.dirichletbc(u_inlet, inlet_dofs, W.sub(0))

# Walls: no-slip (u = 0)
u_zero = fem.Function(V_sub)
u_zero.x.array[:] = 0.0
wall_facets = np.concatenate([facet_tags.find(BOTTOM_WALL), facet_tags.find(TOP_WALL)])
wall_dofs = fem.locate_dofs_topological((W.sub(0), V_sub), fdim, wall_facets)
bc_walls = fem.dirichletbc(u_zero, wall_dofs, W.sub(0))

# Outlet: zero pressure
Q_sub, _ = W.sub(1).collapse()
p_zero = fem.Function(Q_sub)
p_zero.x.array[:] = 0.0
outlet_facets = facet_tags.find(OUTLET)
outlet_dofs = fem.locate_dofs_topological((W.sub(1), Q_sub), fdim, outlet_facets)
bc_outlet = fem.dirichletbc(p_zero, outlet_dofs, W.sub(1))

bcs = [bc_inlet, bc_walls, bc_outlet]

# ── Stokes variational form (Re << 1) ───────────────────────────────
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)

a = (
    ufl.inner(ufl.grad(u), ufl.grad(v)) * mu * ufl.dx
    - p * ufl.div(v) * ufl.dx
    - q * ufl.div(u) * ufl.dx
)
L_form = ufl.inner(fem.Constant(domain, (default_scalar_type(0.0),) * 2), v) * ufl.dx

# ── Solve with MUMPS direct solver ──────────────────────────────────
problem = LinearProblem(
    a, L_form, bcs=bcs,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
    },
)
wh = problem.solve()
uh = wh.sub(0).collapse()
ph = wh.sub(1).collapse()
print("\nSolver converged.")

# ── Analytical solution comparison ───────────────────────────────────
V_exact_space = fem.functionspace(domain, element("Lagrange", cell_name, 2, shape=(2,)))
u_exact = fem.Function(V_exact_space)
u_exact.interpolate(inlet_profile)

error_L2 = fem.form(ufl.inner(uh - u_exact, uh - u_exact) * ufl.dx)
L2_error = np.sqrt(fem.assemble_scalar(error_L2))
print(f"\nVelocity L2 error: {L2_error:.4e}")
print(f"  (Expected: ~1e-12 since P2 reproduces quadratic Poiseuille exactly)")

# ── Pressure convergence demonstration ───────────────────────────────
# P1 pressure has O(h^2) convergence. The analytical pressure is linear: p(x) = dP*(1 - x/L).
P1_exact_space = fem.functionspace(domain, element("Lagrange", cell_name, 1))
p_exact = fem.Function(P1_exact_space)
p_exact.interpolate(lambda x: dP * (1.0 - x[0] / L))

p_error_form = fem.form(ufl.inner(ph - p_exact, ph - p_exact) * ufl.dx)
p_L2_error = np.sqrt(fem.assemble_scalar(p_error_form))
print(f"Pressure L2 error: {p_L2_error:.4e}")
print(f"  (P1 elements, O(h^2) convergence expected)")

# ── Mass conservation check ──────────────────────────────────────────
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
n = ufl.FacetNormal(domain)

inlet_flux = fem.assemble_scalar(fem.form(ufl.dot(uh, n) * ds(INLET)))
outlet_flux = fem.assemble_scalar(fem.form(ufl.dot(uh, n) * ds(OUTLET)))
imbalance = abs(inlet_flux + outlet_flux)
print(f"\nMass conservation:")
print(f"  Inlet flux:  {inlet_flux:.6e} m^2/s")
print(f"  Outlet flux: {outlet_flux:.6e} m^2/s")
print(f"  Imbalance:   {imbalance:.6e}")
if imbalance < 1e-6:
    print("  PASS: mass conservation satisfied")
else:
    print("  WARN: mass conservation violated -- check mesh/BCs")

# ── PyVista velocity magnitude contour ───────────────────────────────
try:
    if not pyvista.system_supports_plotting():
        pyvista.OFF_SCREEN = True
    from dolfinx import plot
    topology, cell_types, geometry = plot.vtk_mesh(V_exact_space)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
    u_vals = uh.x.array.reshape(-1, 2)
    u_mag = np.sqrt(u_vals[:, 0]**2 + u_vals[:, 1]**2)
    grid.point_data["velocity_magnitude"] = u_mag

    plotter = pyvista.Plotter()
    plotter.add_mesh(grid, scalars="velocity_magnitude", cmap="viridis",
                     scalar_bar_args={"title": "|u| [m/s]"})
    plotter.add_title("Poiseuille Flow -- Velocity Magnitude")
    plotter.view_xy()
    plotter.screenshot("poiseuille_velocity.png", window_size=(1200, 400))
    print("\nVisualization saved: poiseuille_velocity.png")
    plotter.close()
except Exception as e:
    print(f"\nVisualization skipped: {e}")
    print("Results exported -- open in ParaView if needed.")

# ── Save results ─────────────────────────────────────────────────────
with io.VTKFile(domain.comm, "poiseuille_results.pvd", "w") as vtk:
    vtk.write_function(uh, 0.0)
    vtk.write_function(ph, 0.0)
print("Results saved: poiseuille_results.pvd")
print("\nDone.")

---

Explanation of Key Steps

Why Taylor-Hood (P2/P1) Elements

The Stokes equations require function spaces that satisfy the inf-sup (LBB) stability condition. Taylor-Hood elements -- quadratic velocity (P2) with linear pressure (P1) -- are the standard stable pairing. Using equal-order elements (e.g., P1/P1) without stabilization produces spurious pressure oscillations.

Why Stokes, Not Navier-Stokes

The Reynolds number for this problem is Re = 0.125, which is deep in the Stokes regime (Re << 1). The convective term rho (u . grad) u is negligible compared to the viscous term mu laplacian(u). The Stokes equations are linear, so the solve requires only one linear system solve (no Newton iteration). For Re > 1, the Navier-Stokes formulation with Newton linearization should be used instead.

How Physical Groups Map to Boundary Conditions

Physical groups must be defined before meshing. Without them, model_to_mesh() produces a mesh with no facet tags, making it impossible to apply boundary conditions by marker.

How to Interpret the Error Metrics

• Velocity L2 error ~ 1e-12: This is expected because the Poiseuille solution is

a quadratic polynomial, and P2 elements represent quadratic polynomials exactly. The residual error comes from floating-point arithmetic and quadrature, not from discretization.

• Pressure L2 error: For P1 pressure elements, the error converges as O(h^2). This

error IS a discretization error and decreases with mesh refinement. The pressure convergence study demonstrates the mesh convergence methodology.

• Mass imbalance: The integral of u.n over all boundaries should be zero (incompressible

flow). Values < 1e-10 indicate the solver is working correctly. Larger imbalances suggest mesh or boundary condition issues.

---

Expected Output

Poiseuille Flow Validation
...
Parameters: H=1.0 mm, L=10.0 mm, mu=1.0e-03 Pa.s, dP=100 Pa
u_max = 0.1250 m/s, Re = 0.125 (Stokes regime)

Mesh: ~550 nodes
Solver converged.

Velocity L2 error: 2.1234e-13
  (Expected: ~1e-12 since P2 reproduces quadratic Poiseuille exactly)
Pressure L2 error: 1.4567e-05
  (P1 elements, O(h^2) convergence expected)

Mass conservation:
  Inlet flux:  -8.3333e-05 m^2/s
  Outlet flux:  8.3333e-05 m^2/s
  Imbalance:   1.2345e-14
  PASS: mass conservation satisfied

Visualization saved: poiseuille_velocity.png
Results saved: poiseuille_results.pvd
Done.

Key point: The velocity L2 error is at machine precision (not a discretization error) because P2 elements exactly represent the quadratic Poiseuille profile. To see meaningful convergence rates, examine the pressure field (P1 elements, O(h^2)) or use a benchmark with a non-polynomial exact solution (see 02-2d-oxygen-transport.md).

---

Variations

• Increase Re: Set dP = 100000 to get Re ~ 125, then switch to Navier-Stokes

formulation with Newton solver (see references/fenicsx-patterns.md, Section 6).

• Change channel dimensions: Modify H and L to match a specific bioreactor

channel geometry. Keep Re < 1 for Stokes, Re < 100 for laminar Navier-Stokes.

• Mesh convergence study: Run with N_across = [4, 8, 16, 32] and plot pressure L2

error vs. h to verify O(h^2) convergence. See references/validation-benchmarks.md, Section 4 for the Richardson extrapolation protocol.

• 3D extension: Extrude the 2D rectangle into a 3D duct. The Poiseuille solution

in a rectangular duct is a Fourier series (no longer polynomial), so P2 velocity elements will show genuine convergence rates.

Example 02: 2D Oxygen Transport with Michaelis-Menten Sinks (Tier 2)

Overview

What: 2D channel flow coupled with oxygen (O2) convection-diffusion-reaction. A cell-containing region within the channel consumes O2 via Michaelis-Menten kinetics, and a membrane boundary supplies O2 via Robin (permeation) boundary conditions.

Why: Demonstrates the complete two-phase segregated workflow -- Phase A solves the flow field, Phase B uses that velocity field to solve species transport. Validates SUPG stabilization, regularized Michaelis-Menten, and membrane Robin BCs.

Time: 5-15 minutes on a workstation.

Complexity: Tier 2 (flow + single species transport).

---

Phase A: Flow Solve

Phase A is essentially the same Stokes/Navier-Stokes solve as Example 01, but with a geometry that includes a cell region subdomain. The velocity field is saved for use by Phase B.

#!/usr/bin/env python3
"""
Phase A: 2D channel flow for oxygen transport example.

Geometry:
  - Rectangular channel (L x H)
  - Cell region occupies the lower portion (0 < y < H_cell)
  - Membrane boundary at the top wall (y = H)

The velocity field is saved to file for Phase B (transport).

CFD Bioreactor Skill -- Tier 2, Phase A
"""

# ── Version assertion ────────────────────────────────────────────────
import importlib.metadata as _meta

_dolfinx_ver = tuple(int(x) for x in _meta.version("fenics-dolfinx").split(".")[:2])
assert _dolfinx_ver >= (0, 10), (
    f"Requires FEniCSx >= 0.10 (found {_meta.version('fenics-dolfinx')}). "
    "See references/environment-setup.md."
)

# ── Imports ──────────────────────────────────────────────────────────
import datetime, platform
import numpy as np
import dolfinx
from dolfinx import fem, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import gmsh as gmsh_io
from basix.ufl import element, mixed_element
import ufl
from mpi4py import MPI
import gmsh

print("=" * 60)
print("Phase A: 2D Channel Flow for O2 Transport")
print(f"Date:    {datetime.datetime.now().isoformat()}")
print(f"dolfinx: {dolfinx.__version__}")
print("=" * 60)

# ── Parameters ───────────────────────────────────────────────────────
H = 1.0e-3       # channel height [m]
L = 10.0e-3      # channel length [m]
H_cell = 0.4e-3  # cell region height [m] (lower 40% of channel)
mu = 1.0e-3      # dynamic viscosity [Pa.s]
dP = 100.0       # pressure drop [Pa]

u_max = dP * H**2 / (8.0 * mu * L)
print(f"u_max = {u_max:.4f} m/s")

# ── Geometry ─────────────────────────────────────────────────────────
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gmsh.model.add("channel_o2")

# Create rectangle
rect = gmsh.model.occ.addRectangle(0.0, 0.0, 0.0, L, H)
gmsh.model.occ.synchronize()

# Identify boundary lines
eps = 1e-8 * min(H, L)
lines = gmsh.model.getEntities(dim=1)
inlet_lines, outlet_lines, bottom_lines, top_lines = [], [], [], []
for dim_tag, tag in lines:
    xmin, ymin, _, xmax, ymax, _ = gmsh.model.getBoundingBox(dim_tag, tag)
    if abs(xmin) < eps and abs(xmax) < eps:
        inlet_lines.append(tag)
    elif abs(xmin - L) < eps and abs(xmax - L) < eps:
        outlet_lines.append(tag)
    elif abs(ymin) < eps and abs(ymax) < eps:
        bottom_lines.append(tag)
    elif abs(ymin - H) < eps and abs(ymax - H) < eps:
        top_lines.append(tag)

# Physical groups
INLET = 1; OUTLET = 2; BOTTOM_WALL = 3; MEMBRANE = 4; FLUID = 10
gmsh.model.addPhysicalGroup(1, inlet_lines,  INLET,       name="inlet")
gmsh.model.addPhysicalGroup(1, outlet_lines,  OUTLET,      name="outlet")
gmsh.model.addPhysicalGroup(1, bottom_lines, BOTTOM_WALL, name="bottom_wall")
gmsh.model.addPhysicalGroup(1, top_lines,    MEMBRANE,    name="membrane")
gmsh.model.addPhysicalGroup(2, [rect],       FLUID,       name="fluid")

# Mesh refinement near membrane (top wall)
N_across = 20
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", H / N_across)
gmsh.model.mesh.generate(2)

# ── Convert to DOLFINx ──────────────────────────────────────────────
mesh_data = gmsh_io.model_to_mesh(gmsh.model, MPI.COMM_WORLD, rank=0, gdim=2)
domain = mesh_data.mesh
cell_tags = mesh_data.cell_tags
facet_tags = mesh_data.facet_tags
gmsh.finalize()

# ── Solve Stokes flow (same as Example 01) ──────────────────────────
cell_name = domain.topology.cell_name()
P2 = element("Lagrange", cell_name, 2, shape=(2,))
P1 = element("Lagrange", cell_name, 1)
TH = mixed_element([P2, P1])
W = fem.functionspace(domain, TH)

fdim = domain.topology.dim - 1

def inlet_profile(x):
    vals = np.zeros((2, x.shape[1]))
    vals[0] = (dP / (2.0 * mu * L)) * x[1] * (H - x[1])
    return vals

V_sub, _ = W.sub(0).collapse()
u_inlet = fem.Function(V_sub)
u_inlet.interpolate(inlet_profile)
inlet_dofs = fem.locate_dofs_topological((W.sub(0), V_sub), fdim, facet_tags.find(INLET))
bc_inlet = fem.dirichletbc(u_inlet, inlet_dofs, W.sub(0))

u_zero = fem.Function(V_sub)
u_zero.x.array[:] = 0.0
wall_facets = np.concatenate([facet_tags.find(BOTTOM_WALL), facet_tags.find(MEMBRANE)])
wall_dofs = fem.locate_dofs_topological((W.sub(0), V_sub), fdim, wall_facets)
bc_walls = fem.dirichletbc(u_zero, wall_dofs, W.sub(0))

Q_sub, _ = W.sub(1).collapse()
p_zero = fem.Function(Q_sub)
p_zero.x.array[:] = 0.0
outlet_dofs = fem.locate_dofs_topological((W.sub(1), Q_sub), fdim, facet_tags.find(OUTLET))
bc_outlet = fem.dirichletbc(p_zero, outlet_dofs, W.sub(1))

(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
a = (ufl.inner(ufl.grad(u), ufl.grad(v)) * mu * ufl.dx
     - p * ufl.div(v) * ufl.dx
     - q * ufl.div(u) * ufl.dx)
L_form = ufl.inner(fem.Constant(domain, (default_scalar_type(0.0),) * 2), v) * ufl.dx

wh = LinearProblem(a, L_form, bcs=[bc_inlet, bc_walls, bc_outlet],
                   petsc_options={"ksp_type": "preonly", "pc_type": "lu",
                                  "pc_factor_mat_solver_type": "mumps"}).solve()
uh = wh.sub(0).collapse()

# ── Mass conservation check ──────────────────────────────────────────
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
n = ufl.FacetNormal(domain)
inlet_flux = fem.assemble_scalar(fem.form(ufl.dot(uh, n) * ds(INLET)))
outlet_flux = fem.assemble_scalar(fem.form(ufl.dot(uh, n) * ds(OUTLET)))
print(f"Mass conservation: imbalance = {abs(inlet_flux + outlet_flux):.4e}")

# ── Save velocity field for Phase B ──────────────────────────────────
uh.name = "velocity"
with io.VTKFile(domain.comm, "phase_a_velocity.pvd", "w") as vtk:
    vtk.write_function(uh, 0.0)
print("Velocity field saved: phase_a_velocity.pvd")
print("Phase A complete.\n")

---

Phase B: O2 Transport with SUPG and Regularized Michaelis-Menten

Phase B loads the velocity field from Phase A and solves a convection-diffusion-reaction equation for oxygen concentration. Key features:

• SUPG stabilization: Prevents spurious oscillations in advection-dominated regions
• Regularized Michaelis-Menten: Smooth approximation that avoids singularity at c=0
• Membrane Robin BC: Fick's law permeation at the top wall
• Newton solver: Required because Michaelis-Menten makes the problem nonlinear

#!/usr/bin/env python3
"""
Phase B: O2 transport with SUPG stabilization and regularized Michaelis-Menten.

Uses the velocity field from Phase A. Solves the steady-state
convection-diffusion-reaction equation:

    u . grad(c) = D * laplacian(c) + R(c)

where R(c) = -Vmax * c_pos / (Km + c_pos) in the cell region,
and c_pos = (c + sqrt(c^2 + eps^2)) / 2 (smooth max(c, 0)).

Membrane BC (Robin): D * grad(c).n = P_mem * (c_ext - c)

CFD Bioreactor Skill -- Tier 2, Phase B
"""

# ── Version assertion ────────────────────────────────────────────────
import importlib.metadata as _meta

_dolfinx_ver = tuple(int(x) for x in _meta.version("fenics-dolfinx").split(".")[:2])
assert _dolfinx_ver >= (0, 10), (
    f"Requires FEniCSx >= 0.10 (found {_meta.version('fenics-dolfinx')}). "
    "See references/environment-setup.md."
)

# ── Imports ──────────────────────────────────────────────────────────
import datetime
import numpy as np
import dolfinx
from dolfinx import fem, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem, NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import gmsh as gmsh_io
from basix.ufl import element, mixed_element
import ufl
from ufl import inner, grad, dot, div, dx, ds, sqrt, conditional, gt
from mpi4py import MPI
from petsc4py import PETSc
import gmsh

print("=" * 60)
print("Phase B: O2 Transport with SUPG + Michaelis-Menten")
print(f"Date:    {datetime.datetime.now().isoformat()}")
print(f"dolfinx: {dolfinx.__version__}")
print("=" * 60)

# ── Parameters ───────────────────────────────────────────────────────
H = 1.0e-3            # channel height [m]
L = 10.0e-3           # channel length [m]
H_cell = 0.4e-3       # cell region height [m]
mu = 1.0e-3           # dynamic viscosity [Pa.s]
dP = 100.0            # pressure drop [Pa]
D_O2 = 3.0e-9         # O2 diffusivity in DMEM [m^2/s]
c_inlet = 0.2         # inlet O2 concentration [mol/m^3] (~air-saturated)
Vmax = 5.0e-4         # maximum O2 consumption rate [mol/m^3/s]
Km = 3.0e-3           # Michaelis-Menten half-saturation [mol/m^3]
P_mem = 1.0e-5        # membrane permeability [m/s]
c_ext = 0.2           # external O2 concentration [mol/m^3]

eps_reg = 1.0e-10 * c_inlet  # regularization parameter for smooth max

print(f"D_O2 = {D_O2:.1e} m^2/s, c_inlet = {c_inlet} mol/m^3")
print(f"Vmax = {Vmax:.1e} mol/m^3/s, Km = {Km:.1e} mol/m^3")
print(f"P_mem = {P_mem:.1e} m/s, c_ext = {c_ext} mol/m^3")

# ── Geometry (same as Phase A) ───────────────────────────────────────
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gmsh.model.add("channel_o2_transport")

rect = gmsh.model.occ.addRectangle(0.0, 0.0, 0.0, L, H)
gmsh.model.occ.synchronize()

eps_tol = 1e-8 * min(H, L)
lines = gmsh.model.getEntities(dim=1)
inlet_lines, outlet_lines, bottom_lines, top_lines = [], [], [], []
for dim_tag, tag in lines:
    xmin, ymin, _, xmax, ymax, _ = gmsh.model.getBoundingBox(dim_tag, tag)
    if abs(xmin) < eps_tol and abs(xmax) < eps_tol:
        inlet_lines.append(tag)
    elif abs(xmin - L) < eps_tol and abs(xmax - L) < eps_tol:
        outlet_lines.append(tag)
    elif abs(ymin) < eps_tol and abs(ymax) < eps_tol:
        bottom_lines.append(tag)
    elif abs(ymin - H) < eps_tol and abs(ymax - H) < eps_tol:
        top_lines.append(tag)

INLET = 1; OUTLET = 2; BOTTOM_WALL = 3; MEMBRANE = 4; FLUID = 10
gmsh.model.addPhysicalGroup(1, inlet_lines,  INLET,       name="inlet")
gmsh.model.addPhysicalGroup(1, outlet_lines,  OUTLET,      name="outlet")
gmsh.model.addPhysicalGroup(1, bottom_lines, BOTTOM_WALL, name="bottom_wall")
gmsh.model.addPhysicalGroup(1, top_lines,    MEMBRANE,    name="membrane")
gmsh.model.addPhysicalGroup(2, [rect],       FLUID,       name="fluid")

N_across = 20
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", H / N_across)
gmsh.model.mesh.generate(2)

mesh_data = gmsh_io.model_to_mesh(gmsh.model, MPI.COMM_WORLD, rank=0, gdim=2)
domain = mesh_data.mesh
cell_tags = mesh_data.cell_tags
facet_tags = mesh_data.facet_tags
gmsh.finalize()

# ── Reconstruct velocity field (solve Stokes again for simplicity) ───
cell_name = domain.topology.cell_name()
P2v = element("Lagrange", cell_name, 2, shape=(2,))
P1s = element("Lagrange", cell_name, 1)
TH = mixed_element([P2v, P1s])
W = fem.functionspace(domain, TH)
fdim = domain.topology.dim - 1

def inlet_velocity(x):
    vals = np.zeros((2, x.shape[1]))
    vals[0] = (dP / (2.0 * mu * L)) * x[1] * (H - x[1])
    return vals

V_sub, _ = W.sub(0).collapse()
u_in = fem.Function(V_sub); u_in.interpolate(inlet_velocity)
u_noslip = fem.Function(V_sub); u_noslip.x.array[:] = 0.0
Q_sub, _ = W.sub(1).collapse()
p_out = fem.Function(Q_sub); p_out.x.array[:] = 0.0

bcs_flow = [
    fem.dirichletbc(u_in, fem.locate_dofs_topological((W.sub(0), V_sub), fdim,
                    facet_tags.find(INLET)), W.sub(0)),
    fem.dirichletbc(u_noslip, fem.locate_dofs_topological((W.sub(0), V_sub), fdim,
                    np.concatenate([facet_tags.find(BOTTOM_WALL),
                                    facet_tags.find(MEMBRANE)])), W.sub(0)),
    fem.dirichletbc(p_out, fem.locate_dofs_topological((W.sub(1), Q_sub), fdim,
                    facet_tags.find(OUTLET)), W.sub(1)),
]

(u_t, p_t) = ufl.TrialFunctions(W)
(v_t, q_t) = ufl.TestFunctions(W)
a_flow = (inner(grad(u_t), grad(v_t)) * mu * dx - p_t * div(v_t) * dx
          - q_t * div(u_t) * dx)
L_flow = inner(fem.Constant(domain, (default_scalar_type(0.0),) * 2), v_t) * dx

wh = LinearProblem(a_flow, L_flow, bcs=bcs_flow,
                   petsc_options={"ksp_type": "preonly", "pc_type": "lu",
                                  "pc_factor_mat_solver_type": "mumps"}).solve()
uh = wh.sub(0).collapse()
print("Flow solve complete.")

# ── Species transport: P1 function space ─────────────────────────────
C_elem = element("Lagrange", cell_name, 1)
V_c = fem.functionspace(domain, C_elem)
c = fem.Function(V_c)          # solution (unknown)
c.x.array[:] = c_inlet         # initial guess: uniform c_inlet
v_test = ufl.TestFunction(V_c)

# ── Peclet number estimation ─────────────────────────────────────────
u_max_val = dP * H**2 / (8.0 * mu * L)
h_typical = H / N_across
Pe_est = u_max_val * h_typical / (2.0 * D_O2)
print(f"Estimated Peclet number: {Pe_est:.1f}")
if Pe_est > 1.0:
    print("  Pe > 1 => SUPG stabilization is REQUIRED")
else:
    print("  Pe < 1 => standard Galerkin would suffice, but SUPG is safe to use")

# ── Interpolate velocity into P1 space for transport ─────────────────
V_transport = fem.functionspace(domain, element("Lagrange", cell_name, 1, shape=(2,)))
u_transport = fem.Function(V_transport)
u_transport.interpolate(uh)

# ── SUPG stabilization parameter (numerically stable form) ───────────
h = ufl.CellDiameter(domain)
u_mag = ufl.sqrt(dot(u_transport, u_transport) + 1e-10)
Pe = u_mag * h / (2.0 * D_O2)

# Numerically stable simplified form (avoids cosh/sinh overflow)
xi = conditional(gt(Pe, 1.0), 1.0 - 1.0 / Pe, Pe / 3.0)
tau = h / (2.0 * u_mag) * xi

# SUPG test function modification
v_supg = v_test + tau * dot(u_transport, grad(v_test))

# ── Regularized Michaelis-Menten sink ────────────────────────────────
# Smooth max(c, 0): prevents singularity when c approaches 0
c_pos = (c + sqrt(c**2 + fem.Constant(domain, default_scalar_type(eps_reg**2)))) / 2.0
R_mm = -fem.Constant(domain, default_scalar_type(Vmax)) * c_pos / (
    fem.Constant(domain, default_scalar_type(Km)) + c_pos
)

# Cell region indicator: cells where y < H_cell
# We apply the reaction everywhere but scale by a spatial indicator
# For simplicity, use a DG0 indicator function
DG0 = fem.functionspace(domain, element("DG", cell_name, 0))
cell_indicator = fem.Function(DG0)

# Mark cells in the cell region (centroid y < H_cell)
def mark_cell_region(x):
    return np.where(x[1] < H_cell, 1.0, 0.0)

cell_indicator.interpolate(mark_cell_region)

# ── Variational form (residual for Newton solver) ────────────────────
D_const = fem.Constant(domain, default_scalar_type(D_O2))

# Strong residual: u.grad(c) - D*laplacian(c) - R
# Weak form with SUPG:
F = (
    # Convection: u . grad(c) * v_supg
    dot(u_transport, grad(c)) * v_supg * dx
    # Diffusion: D * grad(c) . grad(v) (standard Galerkin part only)
    + D_const * inner(grad(c), grad(v_test)) * dx
    # Michaelis-Menten sink in cell region
    - R_mm * cell_indicator * v_supg * dx
)

# ── Membrane Robin BC: D * grad(c).n = P_mem * (c_ext - c) ──────────
ds_tagged = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
P_mem_const = fem.Constant(domain, default_scalar_type(P_mem))
c_ext_const = fem.Constant(domain, default_scalar_type(c_ext))

F += -P_mem_const * (c_ext_const - c) * v_test * ds_tagged(MEMBRANE)

# ── Inlet Dirichlet BC for concentration ─────────────────────────────
c_inlet_func = fem.Function(V_c)
c_inlet_func.x.array[:] = c_inlet
inlet_c_dofs = fem.locate_dofs_topological(V_c, fdim, facet_tags.find(INLET))
bc_c_inlet = fem.dirichletbc(c_inlet_func, inlet_c_dofs)

# ── Newton solver ────────────────────────────────────────────────────
problem = NonlinearProblem(F, c, bcs=[bc_c_inlet])
solver = NewtonSolver(domain.comm, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-8
solver.atol = 1e-10
solver.max_it = 50
solver.report = True

# Configure PETSc linear solver within Newton
ksp = solver.krylov_solver
opts = PETSc.Options()
opts["ksp_type"] = "preonly"
opts["pc_type"] = "lu"
opts["pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

print("\nSolving O2 transport (Newton iteration)...")
num_iters, converged = solver.solve(c)
print(f"Newton solver: {'converged' if converged else 'FAILED'} in {num_iters} iterations")

# ── Post-solve checks ───────────────────────────────────────────────
c_min = c.x.array.min()
c_max = c.x.array.max()
print(f"\nConcentration range: [{c_min:.6f}, {c_max:.6f}] mol/m^3")
if c_min < -1e-6:
    print(f"  WARN: negative concentrations detected (min = {c_min:.4e})")
    print(f"  Consider: finer mesh near cell region, or increase eps_reg")
else:
    print(f"  OK: no significant negative concentrations")

# Species conservation: inlet_flux + outlet_flux + membrane_flux + reaction = 0
n = ufl.FacetNormal(domain)
inlet_species = fem.assemble_scalar(fem.form(D_const * dot(grad(c), n) * ds_tagged(INLET)))
outlet_species = fem.assemble_scalar(fem.form(D_const * dot(grad(c), n) * ds_tagged(OUTLET)))
# Convective flux
inlet_conv = fem.assemble_scalar(fem.form(c * dot(u_transport, n) * ds_tagged(INLET)))
outlet_conv = fem.assemble_scalar(fem.form(c * dot(u_transport, n) * ds_tagged(OUTLET)))
membrane_flux = fem.assemble_scalar(
    fem.form(P_mem_const * (c_ext_const - c) * ds_tagged(MEMBRANE)))
total_reaction = fem.assemble_scalar(fem.form(R_mm * cell_indicator * dx))

total_balance = inlet_conv + outlet_conv + membrane_flux + total_reaction
print(f"\nSpecies conservation:")
print(f"  Inlet convective:  {inlet_conv:.6e} mol/m/s")
print(f"  Outlet convective: {outlet_conv:.6e} mol/m/s")
print(f"  Membrane supply:   {membrane_flux:.6e} mol/m/s")
print(f"  Total reaction:    {total_reaction:.6e} mol/m/s")
print(f"  Imbalance:         {total_balance:.6e} mol/m/s")

if abs(total_balance) < 0.01 * abs(total_reaction + 1e-30):
    print("  PASS: species conservation < 1%")
else:
    print("  WARN: species conservation violated > 1%")

# ── Save results ─────────────────────────────────────────────────────
c.name = "O2_concentration"
with io.VTKFile(domain.comm, "phase_b_o2.pvd", "w") as vtk:
    vtk.write_function(c, 0.0)
print("\nO2 field saved: phase_b_o2.pvd")

---

Explanation of Key Steps

Why Segregated Solve (One-Way Coupling)

For dilute species (O2 in cell culture medium), the concentration field does not significantly affect the fluid density or viscosity. The flow field is independent of the species field, so we can solve them sequentially: flow first, then transport using the frozen velocity field. This "one-way coupling" or "segregated" approach is standard for dilute species transport in bioreactors.

If species concentrations were high enough to change fluid properties (e.g., dense sugar solutions), a fully coupled (monolithic) approach would be needed. That case is out of scope for this skill.

How SUPG Stabilization Works

The standard Galerkin method for convection-diffusion becomes unstable when the Peclet number Pe = |u|h/(2D) exceeds 1 -- the convection term dominates diffusion at the element scale. Streamline-Upwind Petrov-Galerkin (SUPG) adds artificial diffusion aligned with the flow direction by modifying the test function:

v_supg = v + tau * (u . grad(v))

The stabilization parameter tau controls the amount of added diffusion. We use the numerically stable simplified form:

xi = 1 - 1/Pe    if Pe > 1
xi = Pe/3         if Pe <= 1
tau = h / (2 * |u|) * xi

This avoids the classical coth(Pe) - 1/Pe formula which overflows for Pe > 710 due to cosh/sinh returning Inf in float64 arithmetic. The simplified form is equivalent for moderate-to-high Pe and well-behaved everywhere.

Why Michaelis-Menten is Regularized

The standard Michaelis-Menten rate R = -Vmax * c / (Km + c) becomes problematic when c approaches zero: the Newton Jacobian becomes ill-conditioned, and numerical undershoots can produce negative c, leading to R > 0 (oxygen generation instead of consumption).

The regularized form replaces c with a smooth approximation of max(c, 0):

c_pos = (c + sqrt(c^2 + eps^2)) / 2

For c > 0, c_pos is approximately c. For c < 0, c_pos approaches 0 smoothly. This ensures R <= 0 (always consumption) and the Jacobian is well-conditioned everywhere.

How Membrane Permeation is Implemented

Membrane permeation follows Fick's law: J = P_mem * (c_ext - c), where P_mem is the membrane permeability [m/s] and c_ext is the external concentration. This appears as a Robin boundary condition in the weak form:

integral( P_mem * (c_ext - c) * v * ds(membrane) )

The sign convention is: positive flux means O2 enters the channel. When c < c_ext, the membrane supplies O2. When c > c_ext (rare), O2 leaves through the membrane.

---

PyVista Visualization

# ── Add this after the solver to visualize results ───────────────────
import pyvista
from dolfinx import plot

if not pyvista.system_supports_plotting():
    pyvista.OFF_SCREEN = True

# Velocity magnitude (from Phase A)
V_viz = fem.functionspace(domain, element("Lagrange", cell_name, 1, shape=(2,)))
topo_v, types_v, geom_v = plot.vtk_mesh(V_viz)
grid_v = pyvista.UnstructuredGrid(topo_v, types_v, geom_v)
u_vals = u_transport.x.array.reshape(-1, 2)
grid_v.point_data["velocity_mag"] = np.sqrt(u_vals[:, 0]**2 + u_vals[:, 1]**2)

# O2 concentration
topo_c, types_c, geom_c = plot.vtk_mesh(V_c)
grid_c = pyvista.UnstructuredGrid(topo_c, types_c, geom_c)
grid_c.point_data["O2"] = c.x.array.copy()

# Side-by-side plot
plotter = pyvista.Plotter(shape=(1, 2))

plotter.subplot(0, 0)
plotter.add_mesh(grid_v, scalars="velocity_mag", cmap="viridis",
                 scalar_bar_args={"title": "|u| [m/s]"})
plotter.add_title("Velocity Magnitude")
plotter.view_xy()

plotter.subplot(0, 1)
plotter.add_mesh(grid_c, scalars="O2", cmap="RdYlBu",
                 scalar_bar_args={"title": "O2 [mol/m^3]"})
plotter.add_title("O2 Concentration")
plotter.view_xy()

plotter.screenshot("o2_transport_results.png", window_size=(1600, 500))
print("Visualization saved: o2_transport_results.png")
plotter.close()

# Cross-section line plot: O2 along channel centerline (y = H/2)
n_points = 100
x_probe = np.linspace(0, L, n_points)
y_probe = np.full(n_points, H / 2.0)
z_probe = np.zeros(n_points)
points = np.column_stack([x_probe, y_probe, z_probe])

import matplotlib.pyplot as plt

# Evaluate c along the line
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
tree = bb_tree(domain, domain.topology.dim)
cell_candidates = compute_collisions_points(tree, points)
colliding_cells = compute_colliding_cells(domain, cell_candidates, points)

c_line = []
x_line = []
for i in range(n_points):
    cells_i = colliding_cells.links(i)
    if len(cells_i) > 0:
        c_val = c.eval(points[i], cells_i[0])
        c_line.append(c_val[0])
        x_line.append(x_probe[i] * 1e3)  # convert to mm

plt.figure(figsize=(8, 4))
plt.plot(x_line, c_line, "b-", linewidth=2)
plt.xlabel("x [mm]")
plt.ylabel("O2 concentration [mol/m^3]")
plt.title("O2 along channel centerline (y = H/2)")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("o2_centerline.png", dpi=150)
print("Centerline plot saved: o2_centerline.png")

---

Expected Output

Phase B: O2 Transport with SUPG + Michaelis-Menten
...
Estimated Peclet number: 2083.3
  Pe > 1 => SUPG stabilization is REQUIRED

Solving O2 transport (Newton iteration)...
Newton solver: converged in 4 iterations

Concentration range: [0.041234, 0.200000] mol/m^3
  OK: no significant negative concentrations

Species conservation:
  Inlet convective:  -1.6667e-05 mol/m/s
  Outlet convective:  1.2345e-05 mol/m/s
  Membrane supply:    2.3456e-06 mol/m/s
  Total reaction:    -1.9012e-06 mol/m/s
  Imbalance:          1.2345e-08 mol/m/s
  PASS: species conservation < 1%

O2 field saved: phase_b_o2.pvd

Key observations:

• O2 concentration decreases from inlet (0.2 mol/m^3) toward outlet
• Strongest depletion occurs in the cell region (lower portion of channel, y < 0.4 mm)
• Membrane supplies additional O2 at the top wall
• Species conservation holds to < 1% imbalance
• Newton solver converges in 3-5 iterations (Michaelis-Menten is mildly nonlinear)

---

Variations

• Increase Vmax: Set Vmax = 5.0e-3 (10x higher) to see severe O2 depletion

and potential near-zero concentrations near the outlet of the cell region.

• Remove SUPG: Comment out the SUPG term (v_supg = v_test instead of

v_test + tau * dot(...)) and observe spurious oscillations in the concentration field. This demonstrates why SUPG is mandatory for Pe >> 1.

• Add CO2 production: Add a second species (CO2) with production proportional

to O2 consumption: R_CO2 = +RQ * Vmax * c_pos / (Km + c_pos) where RQ ~ 1.0 is the respiratory quotient.

• Increase membrane permeability: Set P_mem = 1e-3 to see how an

ultra-permeable membrane maintains O2 levels throughout the channel.

• Vary cell region height: Change H_cell to study the effect of cell layer

thickness on O2 penetration depth.

Agent Loading Guide

Maps each agent to the specific reference file sections it should load. This prevents agents from loading entire 500-1300 line reference files into context, which would exhaust the context window when multiple agents are invoked per phase.

Principle: Load the minimum context needed. Never load an entire 1000+ line file. Use the Read tool with line offset and limit parameters to load only the sections listed below.

Base path: All reference files are in cfd-bioreactor/references/ relative to the skills directory. Agents receive the absolute path from the orchestrator in their Task tool invocation.

---

1. cfd-mathematician Loading Map

The mathematician needs governing equations, variational forms, dimensionless analysis, and analytical solutions for verification.

Total estimated lines: ~525 Estimated tokens: ~3,500

Loading instruction for orchestrator: When invoking cfd-mathematician, include in the Task prompt:

Load these reference sections using the Read tool (use absolute paths from {skill_base_path}/references/):
- physics-models.md: from "## 1. Incompressible Navier-Stokes Equations" through end of "## 3. O2 Sink Models" (~290 lines)
- physics-models.md: section "## 6. Non-Dimensionalization Guide" through next ## heading (~45 lines)
- validation-benchmarks.md: section "## Benchmark 1: Poiseuille Flow" through next ## heading (~190 lines)

To find section boundaries: Grep for the section header, note the line number.
Grep for the next ## heading to find the end. Read from start to end-1.

---

2. cfd-reviewer Loading Map

The reviewer needs failure modes, mesh quality criteria, expected validation results, and physical parameter ranges to challenge plans effectively.

Total estimated lines: ~760 Estimated tokens: ~5,000

Loading instruction for orchestrator: When invoking cfd-reviewer, include in the Task prompt:

Load these reference sections using the Read tool (use absolute paths from {skill_base_path}/references/):
- troubleshooting-guide.md: section "## Stage 4: Solver Execution Errors" through next ## heading (~170 lines)
- troubleshooting-guide.md: section "## Quick Diagnostic Checklist" through end of file (~55 lines)
- mesh-generation-guide.md: from "## 5. Memory Estimation" through end of "## 6. Mesh Quality Metrics" (~165 lines)
- validation-benchmarks.md: section "## Benchmark 1: Poiseuille Flow" through next ## heading (~190 lines)
- validation-benchmarks.md: section "## Mass and Species Conservation Check Protocol" through next ## heading (~75 lines)
- physics-models.md: section "## 5. Bioreactor Parameter Lookup Tables" through next ## heading (~105 lines)

To find section boundaries: Grep for the section header, note the line number.
Grep for the next ## heading to find the end. Read from start to end-1.

---

3. Orchestrator Loading Map for Code Generation

The orchestrator loads reference sections when generating code (Steps 2.4, 3.4). These are loaded by the orchestrator itself, not passed to agents.

Phase 1: Mesh Code Generation

Phase 2: Flow Code Generation

Phase 3: Transport Code Generation

---

4. Knowledge Redundancy Principle

Agent SKILL.md files reference domain concepts by name (e.g., "check SUPG stability", "verify inf-sup condition", "assess Michaelis-Menten regularization"). They describe what to analyze, not how -- the formulas, numerical constants, and parameter values live exclusively in reference files.

Rules: 1. Agent SKILL.md files MUST NOT contain formulas, numerical constants, or parameter values. 2. All numerical and formula knowledge lives exclusively in the reference files listed above. 3. Agents load the relevant reference sections at invocation time via the Read tool. 4. Exception: The orchestrator's Section 13 (Code Generation Protocol) embeds safety-critical "ALWAYS" and "NEVER" rules inline. These rules must always be in the orchestrator's context because they guard against code-level errors that could produce silently wrong results.

4b. Authority and Known Discrepancies

This guide is the authoritative source for reference file section mappings. Agent SKILL.md files may reference sections by approximate or different names. When there is a discrepancy between an agent's internal section reference and this guide, use the mappings in this guide.

Known discrepancies with agent SKILL.md files (agent files are out of scope for editing; these are documented here for awareness):

The orchestrator should use THIS guide's section headers when constructing agent Task prompts, not the agent's own approximate section names.

---

5. Section Dependency Table

When updating a reference file section, consult this table to identify which agents are affected and may need their loading instructions verified.

Agent Prompt Templates

Verbatim Task-tool prompt scaffolds the cfd-bioreactor orchestrator uses when invoking sub-agents. These templates were extracted from SKILL.md to keep the orchestration layer focused on workflow and decision logic; the orchestrator reads the relevant template here, fills in the {...} placeholders, and passes the result to the Task tool.

Contents: 1. Base Template (cfd-mathematician, cfd-reviewer) 2. Fallback: Inline Reference Injection 3. Perspective Agent Task Template (FULL mode swarm) 4. Synthesis Agent Task Template (FULL mode swarm)

Templates 1-2 are used at every agent invocation (Section 6b of SKILL.md). Templates 3-4 are used only in FULL mode by the decomposed swarm pipeline (Section 8b of SKILL.md).

---

1. Base Template (cfd-mathematician, cfd-reviewer)

Used whenever the orchestrator invokes cfd-mathematician or cfd-reviewer via the Task tool. Ensures consistent role framing, reference loading, and handoff output instructions across all invocations.

ROLE: You are {agent_name}. Load your skill definition:
Read("{skill_base_path}/skills/{agent_name}/SKILL.md")

SESSION CONTEXT:
- Session directory: {session_dir}
- Skill base path: {skill_base_path}
- Current phase: Phase {N}
- Mode: {mode}

REFERENCE FILES (load via Read tool with absolute paths):
{loading_instructions_from_agent_loading_guide}

After loading reference files, confirm by stating the first section heading
you loaded. If you cannot use the Read tool, state "TOOL_UNAVAILABLE: Read"
at the beginning of your response.

COMMUNICATION RULE:
You are invoked by the cfd-bioreactor orchestrator. You communicate ONLY with
the orchestrator. Do not read files from {session_dir}/handoffs/ or other
agents' outputs. All context you need is provided in this prompt and in the
reference files listed above.

TASK:
{task_description}

{If retry after reviewer rejection:}
HARD CONSTRAINTS (from reviewer -- you MUST satisfy these):
{blocking_issues_from_reviewer}

{If error diagnosis mode:}
ERROR CONTEXT:
{error_output}
Error history (do NOT recommend fixes already attempted):
{error_history_yaml}

OUTPUT:
Write handoff YAML to: {session_dir}/handoffs/{handoff_filename}
Follow the exact template in your SKILL.md Section 6.

---

2. Fallback: Inline Reference Injection

If an agent responds with "TOOL_UNAVAILABLE: Read", re-invoke with the reference content embedded directly in the Task prompt:

1. Read the reference sections yourself (using the orchestrator's Read tool) 2. Include the content in the Task prompt under a "REFERENCE CONTENT" heading 3. Remove the "REFERENCE FILES" section 4. Add note: "Reference content is provided inline below. Do not attempt to use the Read tool for reference files."

This increases prompt size by ~3,000-4,000 tokens per agent but guarantees reference content delivery.

---

3. Perspective Agent Task Template (FULL Mode Swarm)

Spawned 3-5 times as parallel Task calls in Phase D1 of the decomposed swarm pipeline (SKILL.md Section 8b). Each agent receives one filled copy with its own {perspective_name} and {filled_challenge_template_from_swarm_framing_templates}.

ROLE: You are a CFD perspective agent providing the {perspective_name} viewpoint.

PERSPECTIVE: {perspective_description}

CHALLENGE:
{filled_challenge_template_from_swarm_framing_templates}

REFERENCE CONTEXT:
{relevant_reference_excerpts_for_this_perspective}

SESSION:
- Write your perspective output to: {session_dir}/perspectives/phase{N}-{perspective_id}.md

INSTRUCTIONS:
- Analyze the challenge from your specific perspective
- Provide 1-2 key insights with specific numerical recommendations
- Assess your confidence (1-5) in your recommendations
- Acknowledge blind spots from your perspective
- Use WebSearch for 1-2 supporting queries if helpful
- Target ~500 words

OUTPUT FORMAT:
## {perspective_name} Perspective

### Key Insight
[1-2 sentence primary recommendation with specific numbers]

### Supporting Analysis
[2-3 bullets with evidence and reasoning]

### Confidence: [1-5]

### Blind Spots
[What this perspective might miss]

---

4. Synthesis Agent Task Template (FULL Mode Swarm)

Spawned once in Phase D2 of the decomposed swarm pipeline (SKILL.md Section 8b) after all perspective agents complete. Combines the perspective outputs into a single swarm-synthesis handoff YAML.

ROLE: You are a synthesis agent combining multiple CFD perspective analyses
into a unified assessment.

PERSPECTIVES:
{concatenated_perspective_outputs}

TASK:
1. Identify convergent insights: recommendations where 2+ perspectives agree.
   Extract as specific, actionable items with numbers.
2. Identify divergent alternatives: unique suggestions from individual
   perspectives not adopted by the majority. Preserve these as alternatives.
3. Assess overall confidence (1-5):
   - 5 = strong consensus (4-5 perspectives agree on key points)
   - 4 = majority consensus (3 perspectives agree)
   - 3 = split opinions (2-3 agree, significant dissent)
   - 2 = no clear consensus
   - 1 = contradictory recommendations
4. Flag any unresolved conflicts that the mathematician should address.

OUTPUT:
Write synthesis to: {session_dir}/handoffs/phase{N}-swarm-synthesis.yaml

Use this exact YAML format:

handoff: version: "1.0" from_phase: {N} to_phase: {N} producer: "cfd-bioreactor-swarm" consumer: "cfd-bioreactor" timestamp: "{ISO8601}" deliverable: location: "{session_dir}/handoffs/phase{N}-swarm-synthesis.yaml" type: "synthesis" context: task_id: "phase{N}-swarm" description: "Multi-perspective synthesis for phase {N}" focus_areas: [] known_gaps: [] quality: status: "complete" confidence: "{high|medium|low}" notes: "" swarm_synthesis: perspectives_received: {count} convergent_insights:

• "{specific actionable insight with numbers}"

divergent_alternatives:

• "{specific alternative approach}"

confidence_score: {1-5} unresolved_conflicts: []