pic.flow

Go to the documentation of this file.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
 
"""!
@file pic.flow
@brief A comprehensive conductor script for the PIC-Flow simulation platform.
 
This script acts as the central user interface for running simulations,
managing configurations, and orchestrating the entire end-to-end workflow.
It translates user-friendly YAML files into C-solver compatible control files,
supports full multi-block configurations, and provides live log streaming.
It features intelligent, content-based config file discovery and robustly
manages data I/O paths for the post-processor. It also supports Slurm job
generation/submission and parameter sweeps via job arrays.
"""
 
import yaml
import sys
import os
import argparse
import subprocess
import shutil
import glob
import csv
import json
import itertools
import re
import shlex
import copy
import numpy as np
from datetime import datetime
import time
 
try:
    import matplotlib.pyplot as plt
except ImportError:
    plt = None
 
# --- Global Path Definitions ---
SCRIPT_PATH = os.path.dirname(os.path.realpath(__file__))
BIN_DIR = SCRIPT_PATH
PROJECT_ROOT = os.path.dirname(BIN_DIR)
 
# Standardized error codes used for CLI/validation reporting.
ERROR_CODE_CLI_USAGE_INVALID = "CLI_USAGE_INVALID"
ERROR_CODE_CFG_MISSING_SECTION = "CFG_MISSING_SECTION"
ERROR_CODE_CFG_MISSING_KEY = "CFG_MISSING_KEY"
ERROR_CODE_CFG_INVALID_TYPE = "CFG_INVALID_TYPE"
ERROR_CODE_CFG_INVALID_VALUE = "CFG_INVALID_VALUE"
ERROR_CODE_CFG_FILE_NOT_FOUND = "CFG_FILE_NOT_FOUND"
ERROR_CODE_CFG_GRID_PARSE = "CFG_GRID_PARSE"
ERROR_CODE_CFG_INCONSISTENT_COMBO = "CFG_INCONSISTENT_COMBO"
 
_ERROR_HINTS = {
    ERROR_CODE_CLI_USAGE_INVALID: "Run 'pic.flow <command> --help' to see valid argument combinations.",
    ERROR_CODE_CFG_MISSING_SECTION: "Add the missing section using examples/master_template/*.yml as reference.",
    ERROR_CODE_CFG_MISSING_KEY: "Add the missing key in the referenced YAML file.",
    ERROR_CODE_CFG_INVALID_TYPE: "Fix the value type to match the documented schema in docs/pages/14_Config_Contract.md.",
    ERROR_CODE_CFG_INVALID_VALUE: "Adjust the value to a supported range/enum from the config reference pages.",
    ERROR_CODE_CFG_FILE_NOT_FOUND: "Fix the path or create the missing file before running again.",
    ERROR_CODE_CFG_GRID_PARSE: "Validate grid file format and numeric payload (block count, dims, coordinates).",
    ERROR_CODE_CFG_INCONSISTENT_COMBO: "Fix conflicting options/keys so the configuration is internally consistent.",
}
 
 
def _sanitize_error_field(value) -> str:
    """Normalize error fields into a single-line string."""
    if value is None:
        return "-"
    text = str(value).strip()
    if not text:
        return "-"
    return " ".join(text.splitlines())
 
 
def emit_structured_error(code: str, key: str = "-", file_path: str = "-",
                          message: str = "", hint: str = None, stream=None):
    """Emit one standardized error line for tooling and users."""
    if stream is None:
        stream = sys.stderr
    resolved_hint = hint if hint is not None else _ERROR_HINTS.get(code, "-")
    print(
        f"ERROR {_sanitize_error_field(code)} | "
        f"key={_sanitize_error_field(key)} | "
        f"file={_sanitize_error_field(file_path)} | "
        f"message={_sanitize_error_field(message)} | "
        f"hint={_sanitize_error_field(resolved_hint)}",
        file=stream,
    )
 
 
def fail_cli_usage(message: str, hint: str = None):
    """Emit a structured CLI usage error and exit with code 2."""
    emit_structured_error(
        ERROR_CODE_CLI_USAGE_INVALID,
        key="-",
        file_path="-",
        message=message,
        hint=hint or _ERROR_HINTS[ERROR_CODE_CLI_USAGE_INVALID],
    )
    sys.exit(2)
 
 
def _split_error_file_and_message(raw_error: str):
    """Split '<file>: <message>' style validation strings when possible."""
    text = str(raw_error).strip()
    match = re.match(r"^(?P<file>[^:]+):\s*(?P<msg>.+)$", text)
    if not match:
        return "-", text
    file_candidate = match.group("file").strip()
    msg = match.group("msg").strip()
    known_suffixes = (".yml", ".yaml", ".cfg", ".picgrid", ".control", ".run", ".txt")
    if "/" in file_candidate or file_candidate.endswith(known_suffixes):
        return file_candidate, msg
    return "-", text
 
 
def _extract_key_path(message: str) -> str:
    """Best-effort key-path extraction from free-form validation messages."""
    dotted = re.search(r"\b([A-Za-z_][A-Za-z0-9_]*(?:\.[A-Za-z0-9_\[\]-]+)+)\b", message)
    if dotted:
        return dotted.group(1)
 
    bracketed = re.search(r"\b([A-Za-z_][A-Za-z0-9_]*\[[^\]]+\](?:\[[^\]]+\])*)\b", message)
    if bracketed:
        return bracketed.group(1)
 
    quoted = re.findall(r"'([A-Za-z0-9_.\[\]-]+)'", message)
    for token in quoted:
        if "." in token or "[" in token or token.isidentifier():
            return token
    return "-"
 
 
def _classify_error_code(message: str) -> str:
    """Map existing validation/error messages to the standardized code set."""
    msg = message.lower()
    if "missing required section" in msg:
        return ERROR_CODE_CFG_MISSING_SECTION
    if "missing required key" in msg or "missing key" in msg:
        return ERROR_CODE_CFG_MISSING_KEY
    if "not found" in msg or "does not exist" in msg:
        return ERROR_CODE_CFG_FILE_NOT_FOUND
    if "invalid dimensions line" in msg or "invalid coordinate row" in msg or "grid file" in msg:
        return ERROR_CODE_CFG_GRID_PARSE
    if (
        "must both be periodic" in msg
        or "inconsistent periodicity" in msg
        or "mismatch" in msg
        or "requires --" in msg
        or "must be 1 (auto) or exactly" in msg
    ):
        return ERROR_CODE_CFG_INCONSISTENT_COMBO
    if (
        "must be a mapping" in msg
        or "must be a list" in msg
        or "must be a string" in msg
        or "must be a boolean" in msg
        or "must be either" in msg
    ):
        return ERROR_CODE_CFG_INVALID_TYPE
    return ERROR_CODE_CFG_INVALID_VALUE
 
# ==============================================================================
# HELPER FUNCTIONS
# ==============================================================================
 
def read_yaml_file(filepath: str) -> dict:
    """!
    @brief Safely reads a YAML file and returns its content.
    @param[in] filepath Path to the YAML file.
    @return A dictionary containing the parsed YAML content.
    @throws SystemExit if the file is not found or cannot be parsed.
    """
    if not os.path.exists(filepath):
        emit_structured_error(
            ERROR_CODE_CFG_FILE_NOT_FOUND,
            key="-",
            file_path=filepath,
            message="Configuration file not found.",
        )
        sys.exit(1)
    try:
        with open(filepath, 'r') as f:
            return yaml.safe_load(f)
    except yaml.YAMLError as e:
        emit_structured_error(
            ERROR_CODE_CFG_INVALID_VALUE,
            key="-",
            file_path=filepath,
            message=f"YAML parse error: {e}",
            hint="Fix YAML syntax/indentation and retry validation.",
        )
        sys.exit(1)
 
def write_yaml_file(filepath: str, data: dict):
    """Write YAML with stable ordering for generated study artifacts."""
    os.makedirs(os.path.dirname(filepath), exist_ok=True)
    with open(filepath, "w") as f:
        yaml.safe_dump(data, f, sort_keys=False)
 
def write_json_file(filepath: str, payload: dict):
    """Write JSON metadata/manifests with a stable, readable format."""
    os.makedirs(os.path.dirname(filepath), exist_ok=True)
    with open(filepath, "w") as f:
        json.dump(payload, f, indent=2, sort_keys=True)
        f.write("\n")
 
def resolve_path(anchor_file: str, candidate: str) -> str:
    """Resolve a potentially relative path against a source YAML file path."""
    if candidate is None:
        return None
    if os.path.isabs(candidate):
        return os.path.abspath(candidate)
    return os.path.abspath(os.path.join(os.path.dirname(os.path.abspath(anchor_file)), candidate))
 
def absolutize_case_external_paths(case_cfg: dict, case_anchor_path: str):
    """Convert external grid/generator paths in case config to absolute paths."""
    grid_cfg = case_cfg.get("grid", {})
    if not isinstance(grid_cfg, dict):
        return
    mode = grid_cfg.get("mode")
    if mode == "file":
        source_file = grid_cfg.get("source_file")
        if isinstance(source_file, str):
            grid_cfg["source_file"] = resolve_path(case_anchor_path, source_file)
    elif mode == "grid_gen":
        gen = grid_cfg.get("generator", {})
        if isinstance(gen, dict):
            for key in ("script", "config_file"):
                val = gen.get(key)
                if isinstance(val, str):
                    gen[key] = resolve_path(case_anchor_path, val)
 
def get_git_commit() -> str:
    """Best-effort git commit lookup for run/study manifests."""
    try:
        result = subprocess.run(
            ["git", "rev-parse", "HEAD"],
            cwd=PROJECT_ROOT,
            text=True,
            capture_output=True,
            check=False
        )
        if result.returncode == 0:
            return result.stdout.strip()
    except Exception:
        pass
    return None
 
def is_valid_email(email: str) -> bool:
    """Lightweight email validation for scheduler notifications."""
    if not isinstance(email, str):
        return False
    pattern = r"^[^@\s]+@[^@\s]+\.[^@\s]+$"
    return re.match(pattern, email.strip()) is not None
 
def normalize_statistics_task(task_name: str) -> str:
    """!
    @brief Normalizes user-facing statistics task names to C pipeline keywords.
    @param[in] task_name Task name from YAML.
    @return Canonical keyword accepted by C statistics pipeline.
    @throws ValueError if task is unsupported.
    """
    if task_name is None:
        raise ValueError("statistics task cannot be None")
    raw = str(task_name).strip()
    if raw == "ComputeMSD":
        return raw
    normalized = raw.lower().replace("-", "_").replace(" ", "_")
    aliases = {
        "msd": "ComputeMSD",
        "compute_msd": "ComputeMSD",
        "computemsd": "ComputeMSD",
    }
    mapped = aliases.get(normalized)
    if mapped is None:
        raise ValueError(f"Unsupported statistics task '{task_name}'. Currently supported: 'msd'.")
    return mapped
 
def _iter_nonempty_noncomment_lines(file_obj):
    """Yield (lineno, stripped_line) for non-empty, non-comment lines."""
    for lineno, raw in enumerate(file_obj, start=1):
        line = raw.strip()
        if not line or line.startswith("#"):
            continue
        yield lineno, line
 
def validate_and_nondimensionalize_picgrid(source_grid: str, dest_grid: str, L_ref: float, expected_nblk: int = None) -> dict:
    """!
    @brief Validates PICGRID payload and writes a non-dimensionalized copy.
    @details Accepts files with or without leading "PICGRID" token. Output is always
             written in canonical PICGRID format with header and per-block dims.
    @param[in] source_grid Input grid file path.
    @param[in] dest_grid Output grid file path.
    @param[in] L_ref Reference length for non-dimensionalization.
    @param[in] expected_nblk Optional expected block count.
    @return Summary dictionary with nblk, dims, and total_nodes.
    @throws ValueError on malformed grid.
    """
    if L_ref == 0.0:
        raise ValueError("length_ref must be non-zero when processing grid coordinates.")
    if not os.path.isfile(source_grid):
        raise ValueError(f"Grid file not found: {source_grid}")
 
    with open(source_grid, "r") as fin:
        line_iter = _iter_nonempty_noncomment_lines(fin)
        try:
            _, first_token = next(line_iter)
        except StopIteration:
            raise ValueError(f"Grid file '{source_grid}' is empty.")
 
        if first_token == "PICGRID":
            try:
                _, nblk_line = next(line_iter)
            except StopIteration:
                raise ValueError(f"Grid file '{source_grid}' missing block count after PICGRID header.")
        else:
            nblk_line = first_token
 
        try:
            nblk = int(nblk_line)
        except ValueError:
            raise ValueError(f"Invalid block count '{nblk_line}' in grid file '{source_grid}'.")
        if nblk <= 0:
            raise ValueError(f"Grid file '{source_grid}' has non-positive block count: {nblk}.")
        if expected_nblk is not None and nblk != expected_nblk:
            raise ValueError(
                f"Grid file block count mismatch: case expects {expected_nblk}, grid contains {nblk}."
            )
 
        dims = []
        for bi in range(nblk):
            try:
                lineno, dim_line = next(line_iter)
            except StopIteration:
                raise ValueError(f"Grid file '{source_grid}' missing dimensions for block {bi}.")
            parts = dim_line.split()
            if len(parts) != 3:
                raise ValueError(
                    f"Invalid dimensions line at {source_grid}:{lineno}. Expected 3 integers, got: '{dim_line}'."
                )
            try:
                im, jm, km = (int(parts[0]), int(parts[1]), int(parts[2]))
            except ValueError:
                raise ValueError(
                    f"Invalid dimensions line at {source_grid}:{lineno}. Non-integer values: '{dim_line}'."
                )
            if im <= 0 or jm <= 0 or km <= 0:
                raise ValueError(
                    f"Invalid block dimensions at {source_grid}:{lineno}: ({im}, {jm}, {km}). Must be > 0."
                )
            dims.append((im, jm, km))
 
        total_nodes_expected = sum(im * jm * km for (im, jm, km) in dims)
        os.makedirs(os.path.dirname(dest_grid), exist_ok=True)
        with open(dest_grid, "w") as fout:
            fout.write("PICGRID\n")
            fout.write(f"{nblk}\n")
            for (im, jm, km) in dims:
                fout.write(f"{im} {jm} {km}\n")
 
            total_nodes_seen = 0
            for lineno, coord_line in line_iter:
                parts = coord_line.split()
                if len(parts) != 3:
                    raise ValueError(
                        f"Invalid coordinate row at {source_grid}:{lineno}. Expected 3 floats, got: '{coord_line}'."
                    )
                try:
                    x = float(parts[0]) / L_ref
                    y = float(parts[1]) / L_ref
                    z = float(parts[2]) / L_ref
                except ValueError:
                    raise ValueError(
                        f"Invalid coordinate row at {source_grid}:{lineno}. Non-numeric values: '{coord_line}'."
                    )
                total_nodes_seen += 1
                if total_nodes_seen > total_nodes_expected:
                    raise ValueError(
                        f"Grid file '{source_grid}' has more coordinates ({total_nodes_seen}) than expected ({total_nodes_expected})."
                    )
                fout.write(f"{x:.8e} {y:.8e} {z:.8e}\n")
 
        if total_nodes_seen != total_nodes_expected:
            raise ValueError(
                f"Grid file '{source_grid}' has {total_nodes_seen} coordinates, expected {total_nodes_expected} from header."
            )
 
    return {"nblk": nblk, "dims": dims, "total_nodes": total_nodes_expected}
 
def run_grid_generator(case_path: str, run_dir: str, grid_cfg: dict) -> str:
    """!
    @brief Runs scripts/grid.gen to produce a PICGRID file for this run.
    @param[in] case_path Path to case.yml (used for relative path resolution).
    @param[in] run_dir Run directory path.
    @param[in] grid_cfg The grid config section from case.yml.
    @return Absolute path to generated dimensional PICGRID file.
    @throws ValueError on invalid config or generator failure.
    """
    generator = grid_cfg.get("generator", {})
    if not isinstance(generator, dict):
        raise ValueError("grid.generator must be a mapping when grid.mode is 'grid_gen'.")
 
    case_dir = os.path.dirname(os.path.abspath(case_path))
    gridgen_script = generator.get("script", os.path.join(SCRIPT_PATH, "grid.gen"))
    if not os.path.isabs(gridgen_script):
        gridgen_script = os.path.abspath(os.path.join(case_dir, gridgen_script))
    if not os.path.isfile(gridgen_script):
        raise ValueError(f"grid.gen script not found: {gridgen_script}")
 
    config_file = generator.get("config_file")
    if not config_file:
        raise ValueError("grid.generator.config_file is required when grid.mode is 'grid_gen'.")
    if not os.path.isabs(config_file):
        config_file = os.path.abspath(os.path.join(case_dir, config_file))
    if not os.path.isfile(config_file):
        raise ValueError(f"grid.generator.config_file not found: {config_file}")
 
    output_file = generator.get("output_file", os.path.join("config", "grid.generated.picgrid"))
    if not os.path.isabs(output_file):
        output_file = os.path.abspath(os.path.join(run_dir, output_file))
    os.makedirs(os.path.dirname(output_file), exist_ok=True)
 
    grid_type = generator.get("grid_type")
    cli_args = generator.get("cli_args", [])
    if cli_args is None:
        cli_args = []
    if not isinstance(cli_args, list):
        raise ValueError("grid.generator.cli_args must be a list of CLI tokens.")
 
    cmd = [sys.executable, gridgen_script, "-c", config_file]
    if grid_type:
        cmd.append(str(grid_type))
    cmd.extend([str(token) for token in cli_args])
    cmd.extend(["--output", output_file])
 
    vts_file = generator.get("vts_file")
    if vts_file:
        if not os.path.isabs(vts_file):
            vts_file = os.path.abspath(os.path.join(run_dir, vts_file))
        os.makedirs(os.path.dirname(vts_file), exist_ok=True)
        cmd.extend(["--vts", vts_file])
 
    stats_file = generator.get("stats_file")
    if stats_file:
        if not os.path.isabs(stats_file):
            stats_file = os.path.abspath(os.path.join(run_dir, stats_file))
        os.makedirs(os.path.dirname(stats_file), exist_ok=True)
        cmd.extend(["--stats-file", stats_file])
 
    print(f"[INFO] Grid generator command: {' '.join(cmd)}")
    result = subprocess.run(cmd, cwd=case_dir, text=True, capture_output=True)
    if result.returncode != 0:
        stderr = (result.stderr or "").strip()
        stdout = (result.stdout or "").strip()
        details = stderr if stderr else stdout
        raise ValueError(
            f"grid.gen failed with exit code {result.returncode}. Details:\n{details}"
        )
    if result.stdout:
        print(result.stdout.strip())
    if result.stderr:
        print(result.stderr.strip())
 
    if not os.path.isfile(output_file):
        raise ValueError(f"grid.gen did not produce expected output file: {output_file}")
 
    return output_file
 
BC_FACE_MAP = {
    "-xi": "-Xi",
    "+xi": "+Xi",
    "-eta": "-Eta",
    "+eta": "+Eta",
    "-zeta": "-Zeta",
    "+zeta": "+Zeta",
}
 
BC_TYPE_MAP = {
    "wall": "WALL",
    "symmetry": "SYMMETRY",
    "inlet": "INLET",
    "outlet": "OUTLET",
    "periodic": "PERIODIC",
}
 
BC_HANDLER_SPECS = {
    # Only handlers that are implemented end-to-end in current C path are allowed.
    "noslip": {
        "types": {"WALL"},
        "required_params": set(),
        "optional_params": set(),
    },
    "constant_velocity": {
        "types": {"INLET"},
        "required_params": {"vx", "vy", "vz"},
        "optional_params": set(),
    },
    "conservation": {
        "types": {"OUTLET"},
        "required_params": set(),
        "optional_params": set(),
    },
    "parabolic": {
        "types": {"INLET"},
        "required_params": {"v_max"},
        "optional_params": set(),
    },
    "geometric": {
        "types": {"PERIODIC"},
        "required_params": set(),
        "optional_params": set(),
    },
    "constant_flux": {
        "types": {"PERIODIC"},
        "required_params": {"target_flux"},
        "optional_params": {"apply_trim"},
    },
}
 
_NUMERIC_BC_PARAMS = {"vx", "vy", "vz", "v_max", "target_flux"}
_BOOL_BC_PARAMS = {"apply_trim"}
 
def _to_float(value, field_name: str) -> float:
    """Convert a YAML scalar to float with a clear error message."""
    try:
        return float(value)
    except (TypeError, ValueError):
        raise ValueError(f"'{field_name}' must be numeric (got {value!r}).")
 
def _to_bool(value, field_name: str) -> bool:
    """Convert a YAML scalar/string to bool with a clear error message."""
    if isinstance(value, bool):
        return value
    if isinstance(value, str):
        raw = value.strip().lower()
        if raw in {"true", "1", "yes"}:
            return True
        if raw in {"false", "0", "no"}:
            return False
    raise ValueError(f"'{field_name}' must be boolean (got {value!r}).")
 
def normalize_boundary_conditions_layout(all_blocks_bcs, num_blocks: int):
    """
    Normalize boundary_conditions to list-of-lists form and validate block count.
    """
    if not all_blocks_bcs:
        raise ValueError("The 'boundary_conditions' section in case.yml is empty.")
 
    is_simple_list = isinstance(all_blocks_bcs[0], dict)
    if num_blocks == 1 and is_simple_list:
        all_blocks_bcs = [all_blocks_bcs]
    elif is_simple_list and num_blocks > 1:
        raise ValueError(
            f"case.yml declares {num_blocks} blocks but boundary_conditions is a single face-list. "
            "Use a list-of-lists, one inner list per block."
        )
 
    if len(all_blocks_bcs) != num_blocks:
        raise ValueError(
            f"Mismatch: case.yml declares {num_blocks} block(s) but found {len(all_blocks_bcs)} BC definitions."
        )
    return all_blocks_bcs
 
def validate_and_prepare_boundary_conditions(case_cfg: dict):
    """
    Validate BC entries against currently supported C-side handlers/types and
    return normalized entries ready for bcs.run generation.
    """
    num_blocks = int(case_cfg.get('models', {}).get('domain', {}).get('blocks', 1))
    scales = case_cfg.get('properties', {}).get('scaling', {})
    L_ref = _to_float(scales.get('length_ref'), "properties.scaling.length_ref")
    U_ref = _to_float(scales.get('velocity_ref'), "properties.scaling.velocity_ref")
    if U_ref == 0.0:
        raise ValueError("properties.scaling.velocity_ref must be non-zero for non-dimensionalization.")
    if L_ref == 0.0:
        raise ValueError("properties.scaling.length_ref must be non-zero for non-dimensionalization.")
 
    all_blocks_bcs = normalize_boundary_conditions_layout(case_cfg.get('boundary_conditions', []), num_blocks)
    prepared_blocks = []
 
    expected_faces = {"-Xi", "+Xi", "-Eta", "+Eta", "-Zeta", "+Zeta"}
    axis_pairs = [("-Xi", "+Xi"), ("-Eta", "+Eta"), ("-Zeta", "+Zeta")]
 
    for bi, block_bcs in enumerate(all_blocks_bcs):
        if not isinstance(block_bcs, list):
            raise ValueError(f"boundary_conditions[{bi}] must be a list of face configs.")
 
        prepared_block = []
        seen_faces = {}
 
        for idx, bc in enumerate(block_bcs):
            if not isinstance(bc, dict):
                raise ValueError(f"boundary_conditions[{bi}][{idx}] must be a mapping.")
 
            for req in ("face", "type", "handler"):
                if req not in bc:
                    raise ValueError(f"boundary_conditions[{bi}][{idx}] missing required key '{req}'.")
 
            face_raw = str(bc["face"]).strip()
            face_key = face_raw.lower()
            face = BC_FACE_MAP.get(face_key)
            if face is None:
                raise ValueError(
                    f"Unsupported BC face '{face_raw}' at boundary_conditions[{bi}][{idx}]. "
                    f"Supported: {sorted(expected_faces)}."
                )
            if face in seen_faces:
                raise ValueError(f"Duplicate face '{face}' in boundary_conditions[{bi}] (entries {seen_faces[face]} and {idx}).")
            seen_faces[face] = idx
 
            bc_type_raw = str(bc["type"]).strip()
            bc_type = BC_TYPE_MAP.get(bc_type_raw.lower())
            if bc_type is None:
                raise ValueError(
                    f"Unsupported BC type '{bc_type_raw}' for face {face} in block {bi}. "
                    f"Supported: {sorted(set(BC_TYPE_MAP.values()))}."
                )
 
            handler = str(bc["handler"]).strip().lower()
            handler_spec = BC_HANDLER_SPECS.get(handler)
            if handler_spec is None:
                raise ValueError(
                    f"Unsupported BC handler '{bc['handler']}' for face {face} in block {bi}. "
                    f"Supported now: {sorted(BC_HANDLER_SPECS.keys())}."
                )
            if bc_type not in handler_spec["types"]:
                raise ValueError(
                    f"Invalid BC combination on block {bi}, face {face}: type '{bc_type}' cannot use handler '{handler}'."
                )
 
            params = bc.get("params", {})
            if params is None:
                params = {}
            if not isinstance(params, dict):
                raise ValueError(f"'params' for block {bi}, face {face} must be a mapping.")
 
            # Reject legacy structured keys explicitly.
            if "vector" in params or "velocity" in params:
                raise ValueError(
                    f"Legacy params key ('vector'/'velocity') found on block {bi}, face {face}. "
                    "Use scalar keys 'vx', 'vy', 'vz'."
                )
 
            required = handler_spec["required_params"]
            optional = handler_spec["optional_params"]
            allowed = required | optional
 
            missing = sorted(required - set(params.keys()))
            if missing:
                raise ValueError(
                    f"Missing required params for handler '{handler}' on block {bi}, face {face}: {missing}."
                )
            unknown = sorted(set(params.keys()) - allowed)
            if unknown:
                raise ValueError(
                    f"Unknown params for handler '{handler}' on block {bi}, face {face}: {unknown}. "
                    f"Allowed: {sorted(allowed)}."
                )
 
            converted_params = {}
            for key, value in params.items():
                if key in _NUMERIC_BC_PARAMS:
                    numeric = _to_float(value, f"boundary_conditions[{bi}][{idx}].params.{key}")
                    if key in {"vx", "vy", "vz", "v_max"}:
                        converted_params[key] = numeric / U_ref
                    elif key == "target_flux":
                        converted_params[key] = numeric / (U_ref * (L_ref ** 2))
                elif key in _BOOL_BC_PARAMS:
                    converted_params[key] = _to_bool(value, f"boundary_conditions[{bi}][{idx}].params.{key}")
                else:
                    # Defensive fallback; should not happen due unknown-key gate above.
                    converted_params[key] = value
 
            prepared_block.append({
                "face": face,
                "type": bc_type,
                "handler": handler,
                "params": converted_params,
            })
 
        missing_faces = sorted(expected_faces - set(seen_faces.keys()))
        if missing_faces:
            raise ValueError(
                f"boundary_conditions[{bi}] is incomplete. Missing faces: {missing_faces}. "
                "Provide all six faces explicitly."
            )
 
        # Pairwise periodic consistency checks.
        face_map = {entry["face"]: entry for entry in prepared_block}
        for neg_face, pos_face in axis_pairs:
            neg = face_map[neg_face]
            pos = face_map[pos_face]
            neg_periodic = (neg["type"] == "PERIODIC")
            pos_periodic = (pos["type"] == "PERIODIC")
            if neg_periodic != pos_periodic:
                raise ValueError(
                    f"Inconsistent periodicity in block {bi}: {neg_face} and {pos_face} must both be PERIODIC or neither."
                )
 
            driven_handlers = {"constant_flux"}
            if (neg["handler"] in driven_handlers) or (pos["handler"] in driven_handlers):
                if neg["handler"] != pos["handler"]:
                    raise ValueError(
                        f"In block {bi}, driven periodic handlers on {neg_face}/{pos_face} must match exactly."
                    )
                if not (neg_periodic and pos_periodic):
                    raise ValueError(
                        f"In block {bi}, driven periodic handler '{neg['handler']}' requires PERIODIC type on both faces."
                    )
 
        prepared_blocks.append(prepared_block)
 
    return prepared_blocks
 
def validate_solver_configs(case_cfg: dict, solver_cfg: dict, monitor_cfg: dict,
                            case_path: str, solver_path: str, monitor_path: str):
    """!
    @brief Validates all solver input configs before any work is done.
    @details Checks for required sections, required keys, and physical sanity.
             Exits with a clear error message on the first problem found.
    @param[in] case_cfg    Parsed case YAML dictionary.
    @param[in] solver_cfg  Parsed solver YAML dictionary.
    @param[in] monitor_cfg Parsed monitor YAML dictionary.
    @param[in] case_path   Path to case file (for error messages).
    @param[in] solver_path Path to solver file (for error messages).
    @param[in] monitor_path Path to monitor file (for error messages).
    @throws SystemExit on validation failure.
    """
    errors = []
 
    # --- case.yml: required top-level sections ---
    required_case_sections = ['properties', 'run_control', 'grid', 'models', 'boundary_conditions']
    for section in required_case_sections:
        if section not in case_cfg:
            errors.append(f"  {case_path}: missing required section '{section}'.")
 
    if errors:
        _print_validation_errors(errors)
 
    # --- case.yml: properties sub-keys ---
    props = case_cfg.get('properties', {})
    for group, keys in [('scaling', ['length_ref', 'velocity_ref']),
                        ('fluid', ['density', 'viscosity']),
                        ('initial_conditions', ['u_physical', 'v_physical', 'w_physical'])]:
        sub = props.get(group, {})
        if not sub:
            errors.append(f"  {case_path}: missing 'properties.{group}' section.")
        else:
            for k in keys:
                if k not in sub:
                    errors.append(f"  {case_path}: missing key 'properties.{group}.{k}'.")
            if group == 'initial_conditions' and 'mode' in sub:
                try:
                    normalize_field_init_mode(sub.get('mode'))
                except ValueError as e:
                    errors.append(f"  {case_path}: {e}")
 
    # --- case.yml: run_control sub-keys ---
    rc = case_cfg.get('run_control', {})
    for k in ['start_step', 'total_steps', 'dt_physical']:
        if k not in rc:
            errors.append(f"  {case_path}: missing key 'run_control.{k}'.")
 
    # --- Physical sanity checks ---
    try:
        density = float(props.get('fluid', {}).get('density', 0))
        viscosity = float(props.get('fluid', {}).get('viscosity', 0))
        dt = float(rc.get('dt_physical', 0))
        if density <= 0:
            errors.append(f"  {case_path}: 'properties.fluid.density' must be positive (got {density}).")
        if viscosity < 0:
            errors.append(f"  {case_path}: 'properties.fluid.viscosity' must be non-negative (got {viscosity}).")
        if dt <= 0:
            errors.append(f"  {case_path}: 'run_control.dt_physical' must be positive (got {dt}).")
    except (TypeError, ValueError):
        pass  # Will be caught later during processing
 
    # --- case.yml: grid mode ---
    grid_cfg = case_cfg.get('grid', {})
    grid_mode = grid_cfg.get('mode')
    valid_grid_modes = ['file', 'programmatic_c', 'grid_gen']
    if grid_mode not in valid_grid_modes:
        errors.append(f"  {case_path}: 'grid.mode' must be one of {valid_grid_modes} (got '{grid_mode}').")
    elif grid_mode == 'file':
        source_file = grid_cfg.get('source_file')
        if not source_file:
            errors.append(f"  {case_path}: 'grid.source_file' is required when grid.mode is 'file'.")
        else:
            source_abs = source_file if os.path.isabs(source_file) else os.path.abspath(os.path.join(os.path.dirname(case_path), source_file))
            if not os.path.isfile(source_abs):
                errors.append(f"  {case_path}: grid.source_file does not exist: {source_abs}")
    elif grid_mode == 'programmatic_c':
        grid_settings = grid_cfg.get('programmatic_settings')
        if not grid_settings:
            errors.append(f"  {case_path}: 'grid.programmatic_settings' is required when grid.mode is 'programmatic_c'.")
        elif not isinstance(grid_settings, dict):
            errors.append(f"  {case_path}: 'grid.programmatic_settings' must be a mapping.")
        else:
            for p_key in ('da_processors_x', 'da_processors_y', 'da_processors_z'):
                p_val = grid_settings.get(p_key)
                if isinstance(p_val, (list, tuple)):
                    errors.append(
                        f"  {case_path}: grid.programmatic_settings.{p_key} must be a scalar integer. "
                        "Per-block MPI decomposition is not implemented on the C side; DMDA layout is global."
                    )
                elif p_val is not None and (not isinstance(p_val, int) or p_val <= 0):
                    errors.append(
                        f"  {case_path}: grid.programmatic_settings.{p_key} must be a positive integer when provided (got {p_val})."
                    )
    elif grid_mode == 'grid_gen':
        gen_cfg = grid_cfg.get('generator')
        if not isinstance(gen_cfg, dict):
            errors.append(f"  {case_path}: 'grid.generator' must be a mapping when grid.mode is 'grid_gen'.")
        else:
            config_file = gen_cfg.get('config_file')
            if not config_file:
                errors.append(f"  {case_path}: 'grid.generator.config_file' is required for grid.mode='grid_gen'.")
            else:
                config_abs = config_file if os.path.isabs(config_file) else os.path.abspath(os.path.join(os.path.dirname(case_path), config_file))
                if not os.path.isfile(config_abs):
                    errors.append(f"  {case_path}: grid.generator.config_file does not exist: {config_abs}")
 
            grid_type = gen_cfg.get('grid_type')
            if grid_type is not None and str(grid_type) not in {'cpipe', 'pipe', 'warp'}:
                errors.append(f"  {case_path}: grid.generator.grid_type must be one of ['cpipe','pipe','warp'] (got '{grid_type}').")
 
            cli_args = gen_cfg.get('cli_args', [])
            if cli_args is not None and not isinstance(cli_args, list):
                errors.append(f"  {case_path}: grid.generator.cli_args must be a list of CLI tokens.")
 
    # --- case.yml: boundary_conditions strict validation ---
    try:
        validate_and_prepare_boundary_conditions(case_cfg)
    except ValueError as e:
        errors.append(f"  {case_path}: {e}")
 
    # --- case.yml: particle initialization validation ---
    particles_cfg = case_cfg.get('models', {}).get('physics', {}).get('particles', {})
    if particles_cfg and not isinstance(particles_cfg, dict):
        errors.append(f"  {case_path}: 'models.physics.particles' must be a mapping.")
    elif isinstance(particles_cfg, dict):
        init_mode_raw = particles_cfg.get('init_mode', 'Surface')
        try:
            pinit_code = normalize_particle_init_mode(init_mode_raw)
        except ValueError as e:
            errors.append(f"  {case_path}: {e}")
            pinit_code = None
 
        restart_mode = particles_cfg.get('restart_mode')
        if restart_mode is not None and str(restart_mode).lower() not in {"init", "load"}:
            errors.append(
                f"  {case_path}: models.physics.particles.restart_mode must be 'init' or 'load' (got '{restart_mode}')."
            )
 
        if pinit_code == 2:
            point_cfg = particles_cfg.get('point_source', {})
            if not isinstance(point_cfg, dict):
                errors.append(f"  {case_path}: models.physics.particles.point_source must be a mapping when init_mode is PointSource.")
            else:
                for coord in ('x', 'y', 'z'):
                    if coord not in point_cfg:
                        errors.append(
                            f"  {case_path}: models.physics.particles.point_source.{coord} is required when init_mode is PointSource."
                        )
 
    # --- solver.yml: basic structure ---
    if not isinstance(solver_cfg, dict) or not solver_cfg:
        errors.append(f"  {solver_path}: solver config is empty or not a valid YAML mapping.")
    else:
        strategy_cfg = solver_cfg.get('strategy', {})
        if not isinstance(strategy_cfg, dict):
            errors.append(f"  {solver_path}: 'strategy' must be a mapping.")
        elif 'implicit' in strategy_cfg:
            errors.append(
                f"  {solver_path}: legacy key 'strategy.implicit' is not supported. "
                "Use 'strategy.momentum_solver' with named solver values."
            )
        if isinstance(strategy_cfg, dict) and 'momentum_solver' in strategy_cfg:
            try:
                normalize_momentum_solver_type(strategy_cfg['momentum_solver'])
            except ValueError as e:
                errors.append(f"  {solver_path}: {e}")
 
        op_mode_cfg = solver_cfg.get('operation_mode', {})
        if op_mode_cfg is not None and not isinstance(op_mode_cfg, dict):
            errors.append(f"  {solver_path}: 'operation_mode' must be a mapping when provided.")
        elif isinstance(op_mode_cfg, dict):
            analytical_type = op_mode_cfg.get('analytical_type')
            if analytical_type is not None and not isinstance(analytical_type, str):
                errors.append(f"  {solver_path}: 'operation_mode.analytical_type' must be a string when provided.")
 
        ms_cfg = solver_cfg.get('momentum_solver', {})
        if ms_cfg is not None and not isinstance(ms_cfg, dict):
            errors.append(f"  {solver_path}: 'momentum_solver' must be a mapping when provided.")
        elif isinstance(ms_cfg, dict):
            legacy_flat_keys = {
                'max_pseudo_steps', 'absolute_tol', 'relative_tol', 'step_tol',
                'pseudo_cfl', 'rk4_residual_noise_allowance_factor'
            }
            present_legacy = sorted(legacy_flat_keys.intersection(ms_cfg.keys()))
            if present_legacy:
                errors.append(
                    f"  {solver_path}: legacy flat keys in 'momentum_solver' are not supported: {present_legacy}. "
                    "Use solver-specific sub-blocks (e.g., momentum_solver.dual_time_picard_rk4)."
                )
 
            if 'type' in ms_cfg:
                try:
                    normalize_momentum_solver_type(ms_cfg['type'])
                except ValueError as e:
                    errors.append(f"  {solver_path}: {e}")
 
            allowed_ms_keys = {'type', 'dual_time_picard_rk4'}
            unknown_ms_keys = sorted(set(ms_cfg.keys()) - allowed_ms_keys)
            if unknown_ms_keys:
                errors.append(
                    f"  {solver_path}: unsupported momentum_solver blocks/keys: {unknown_ms_keys}. "
                    "Currently supported: 'dual_time_picard_rk4' (plus optional 'type')."
                )
 
            selected_solver = None
            if isinstance(strategy_cfg, dict) and 'momentum_solver' in strategy_cfg:
                try:
                    selected_solver = normalize_momentum_solver_type(strategy_cfg['momentum_solver'])
                except ValueError:
                    pass
            if selected_solver is None and 'type' in ms_cfg:
                try:
                    selected_solver = normalize_momentum_solver_type(ms_cfg['type'])
                except ValueError:
                    pass
            if selected_solver is None:
                selected_solver = "DUALTIME_PICARD_RK4"
 
            if selected_solver != "DUALTIME_PICARD_RK4" and 'dual_time_picard_rk4' in ms_cfg:
                errors.append(
                    f"  {solver_path}: momentum_solver.dual_time_picard_rk4 is set but selected solver is "
                    f"{selected_solver}. Solver-specific blocks must match the selected solver."
                )
 
            dt_picard_cfg = ms_cfg.get('dual_time_picard_rk4')
            if dt_picard_cfg is not None:
                if not isinstance(dt_picard_cfg, dict):
                    errors.append(f"  {solver_path}: momentum_solver.dual_time_picard_rk4 must be a mapping.")
                else:
                    allowed_dt_keys = {
                        'max_pseudo_steps', 'absolute_tol', 'relative_tol', 'step_tol',
                        'pseudo_cfl', 'rk4_residual_noise_allowance_factor'
                    }
                    unknown_dt_keys = sorted(set(dt_picard_cfg.keys()) - allowed_dt_keys)
                    if unknown_dt_keys:
                        errors.append(
                            f"  {solver_path}: unsupported keys in momentum_solver.dual_time_picard_rk4: {unknown_dt_keys}."
                        )
                    if 'pseudo_cfl' in dt_picard_cfg:
                        pcfl_cfg = dt_picard_cfg['pseudo_cfl']
                        if not isinstance(pcfl_cfg, dict):
                            errors.append(f"  {solver_path}: momentum_solver.dual_time_picard_rk4.pseudo_cfl must be a mapping.")
                        else:
                            allowed_pcfl_keys = {'initial', 'minimum', 'maximum', 'growth_factor', 'reduction_factor'}
                            unknown_pcfl_keys = sorted(set(pcfl_cfg.keys()) - allowed_pcfl_keys)
                            if unknown_pcfl_keys:
                                errors.append(
                                    f"  {solver_path}: unsupported keys in momentum_solver.dual_time_picard_rk4.pseudo_cfl: {unknown_pcfl_keys}."
                                )
 
    # --- monitor.yml: basic structure ---
    if not isinstance(monitor_cfg, dict) or not monitor_cfg:
        errors.append(f"  {monitor_path}: monitor config is empty or not a valid YAML mapping.")
    else:
        io_cfg = monitor_cfg.get('io', {})
        freq = io_cfg.get('data_output_frequency')
        if freq is not None and (not isinstance(freq, int) or freq <= 0):
            errors.append(f"  {monitor_path}: 'io.data_output_frequency' must be a positive integer (got {freq}).")
        solver_mon_cfg = monitor_cfg.get('solver_monitoring')
        if solver_mon_cfg is not None and not isinstance(solver_mon_cfg, dict):
            errors.append(f"  {monitor_path}: 'solver_monitoring' must be a mapping of <flag>: <value>.")
 
    if errors:
        _print_validation_errors(errors)
 
 
def validate_post_config(post_cfg: dict, post_path: str):
    """!
    @brief Validates the post-processing config before running the post-processor.
    @param[in] post_cfg  Parsed post-processing YAML dictionary.
    @param[in] post_path Path to post file (for error messages).
    @throws SystemExit on validation failure.
    """
    errors = []
 
    if not isinstance(post_cfg, dict) or not post_cfg:
        errors.append(f"  {post_path}: post-processing config is empty or not a valid YAML mapping.")
        _print_validation_errors(errors)
 
    # --- run_control ---
    if 'run_control' not in post_cfg:
        errors.append(f"  {post_path}: missing required section 'run_control'.")
 
    # --- io section ---
    io_cfg = post_cfg.get('io', {})
    if not io_cfg:
        errors.append(f"  {post_path}: missing required section 'io'.")
    else:
        for k in ['output_directory', 'output_filename_prefix']:
            if k not in io_cfg:
                errors.append(f"  {post_path}: missing required key 'io.{k}'.")
        input_extensions = io_cfg.get('input_extensions')
        if input_extensions is not None:
            if not isinstance(input_extensions, dict):
                errors.append(f"  {post_path}: 'io.input_extensions' must be a mapping when provided.")
            else:
                for ext_key in ('eulerian', 'particle'):
                    ext_val = input_extensions.get(ext_key)
                    if ext_val is not None and not isinstance(ext_val, str):
                        errors.append(f"  {post_path}: 'io.input_extensions.{ext_key}' must be a string extension.")
 
        averaged_fields = io_cfg.get('eulerian_fields_averaged')
        if averaged_fields is not None and not isinstance(averaged_fields, list):
            errors.append(f"  {post_path}: 'io.eulerian_fields_averaged' must be a list when provided.")
 
    # --- Check eulerian_pipeline entries have 'task' key ---
    for i, entry in enumerate(post_cfg.get('eulerian_pipeline', [])):
        if not isinstance(entry, dict) or 'task' not in entry:
            errors.append(f"  {post_path}: 'eulerian_pipeline[{i}]' is missing the 'task' key. "
                          "Check YAML indentation (each entry needs '- task: ...' with proper spacing).")
 
    # --- Check statistics pipeline entries ---
    stats_cfg = post_cfg.get('statistics_pipeline')
    stats_entries = []
    if stats_cfg is not None:
        if isinstance(stats_cfg, list):
            stats_entries = stats_cfg
        elif isinstance(stats_cfg, dict):
            stats_entries = stats_cfg.get('tasks', [])
            if not isinstance(stats_entries, list):
                errors.append(f"  {post_path}: 'statistics_pipeline.tasks' must be a list.")
            stats_output_prefix = stats_cfg.get('output_prefix')
            if stats_output_prefix is not None and not isinstance(stats_output_prefix, str):
                errors.append(f"  {post_path}: 'statistics_pipeline.output_prefix' must be a string.")
        else:
            errors.append(
                f"  {post_path}: 'statistics_pipeline' must be either a list of tasks or a mapping with a 'tasks' list."
            )
        for i, entry in enumerate(stats_entries):
            if isinstance(entry, str):
                task_name = entry
            elif isinstance(entry, dict) and 'task' in entry:
                task_name = entry.get('task')
            else:
                errors.append(
                    f"  {post_path}: statistics task entry {i} must be either a string or a mapping with key 'task'."
                )
                continue
            try:
                normalize_statistics_task(task_name)
            except ValueError as e:
                errors.append(f"  {post_path}: {e}")
 
    if errors:
        _print_validation_errors(errors)
 
def validate_cluster_config(cluster_cfg: dict, cluster_path: str):
    """Validate Slurm scheduler configuration from cluster.yml."""
    errors = []
    if not isinstance(cluster_cfg, dict) or not cluster_cfg:
        errors.append(f"  {cluster_path}: cluster config is empty or not a valid YAML mapping.")
        _print_validation_errors(errors)
 
    scheduler = cluster_cfg.get("scheduler", {})
    if not isinstance(scheduler, dict):
        errors.append(f"  {cluster_path}: 'scheduler' must be a mapping.")
    else:
        scheduler_type = scheduler.get("type", "slurm")
        if str(scheduler_type).lower() != "slurm":
            errors.append(f"  {cluster_path}: scheduler.type must be 'slurm' in v1 (got '{scheduler_type}').")
 
    resources = cluster_cfg.get("resources", {})
    if not isinstance(resources, dict):
        errors.append(f"  {cluster_path}: 'resources' must be a mapping.")
    else:
        for req in ("account", "nodes", "ntasks_per_node", "mem", "time"):
            if req not in resources:
                errors.append(f"  {cluster_path}: missing required key 'resources.{req}'.")
        for int_key in ("nodes", "ntasks_per_node"):
            if int_key in resources:
                val = resources.get(int_key)
                if not isinstance(val, int) or val <= 0:
                    errors.append(f"  {cluster_path}: resources.{int_key} must be a positive integer (got {val}).")
        for str_key in ("account", "mem", "time", "partition"):
            if str_key in resources and resources.get(str_key) is not None:
                if not isinstance(resources.get(str_key), str):
                    errors.append(f"  {cluster_path}: resources.{str_key} must be a string when provided.")
 
    notifications = cluster_cfg.get("notifications", {})
    if notifications is not None and not isinstance(notifications, dict):
        errors.append(f"  {cluster_path}: 'notifications' must be a mapping when provided.")
    elif isinstance(notifications, dict):
        mail_user = notifications.get("mail_user")
        if mail_user is not None and not is_valid_email(mail_user):
            errors.append(f"  {cluster_path}: notifications.mail_user is not a valid email '{mail_user}'.")
        mail_type = notifications.get("mail_type")
        if mail_type is not None and not isinstance(mail_type, str):
            errors.append(f"  {cluster_path}: notifications.mail_type must be a string when provided.")
 
    execution = cluster_cfg.get("execution", {})
    if execution is not None and not isinstance(execution, dict):
        errors.append(f"  {cluster_path}: 'execution' must be a mapping when provided.")
    elif isinstance(execution, dict):
        module_setup = execution.get("module_setup", [])
        if module_setup is not None and not isinstance(module_setup, list):
            errors.append(f"  {cluster_path}: execution.module_setup must be a list of shell lines.")
        elif isinstance(module_setup, list):
            for i, line in enumerate(module_setup):
                if not isinstance(line, str):
                    errors.append(f"  {cluster_path}: execution.module_setup[{i}] must be a string.")
 
        launcher = execution.get("launcher", "srun")
        if launcher is not None and not isinstance(launcher, str):
            errors.append(f"  {cluster_path}: execution.launcher must be a string when provided.")
        launcher_args = execution.get("launcher_args", [])
        if launcher_args is not None and not isinstance(launcher_args, list):
            errors.append(f"  {cluster_path}: execution.launcher_args must be a list of CLI tokens.")
        elif isinstance(launcher_args, list):
            for i, token in enumerate(launcher_args):
                if not isinstance(token, (str, int, float)):
                    errors.append(f"  {cluster_path}: execution.launcher_args[{i}] must be a scalar CLI token.")
 
        extra_sbatch = execution.get("extra_sbatch")
        if extra_sbatch is not None and not isinstance(extra_sbatch, (dict, list)):
            errors.append(f"  {cluster_path}: execution.extra_sbatch must be a mapping or list when provided.")
 
    if errors:
        _print_validation_errors(errors)
 
def validate_study_config(study_cfg: dict, study_path: str):
    """Validate sweep/study specification from study.yml."""
    errors = []
    if not isinstance(study_cfg, dict) or not study_cfg:
        errors.append(f"  {study_path}: study config is empty or not a valid YAML mapping.")
        _print_validation_errors(errors)
 
    base_cfgs = study_cfg.get("base_configs")
    if not isinstance(base_cfgs, dict):
        errors.append(f"  {study_path}: missing required mapping 'base_configs'.")
    else:
        for req in ("case", "solver", "monitor", "post"):
            path_val = base_cfgs.get(req)
            if not path_val or not isinstance(path_val, str):
                errors.append(f"  {study_path}: base_configs.{req} must be a path string.")
            else:
                resolved = resolve_path(study_path, path_val)
                if not os.path.isfile(resolved):
                    errors.append(f"  {study_path}: base_configs.{req} does not exist: {resolved}")
 
    study_type = study_cfg.get("study_type")
    allowed_types = {"grid_independence", "timestep_independence", "sensitivity"}
    if study_type not in allowed_types:
        errors.append(
            f"  {study_path}: study_type must be one of {sorted(allowed_types)} (got '{study_type}')."
        )
 
    parameters = study_cfg.get("parameters")
    if not isinstance(parameters, dict) or not parameters:
        errors.append(f"  {study_path}: 'parameters' must be a non-empty mapping of key->list.")
    else:
        allowed_roots = {"case", "solver", "monitor", "post"}
        for key, values in parameters.items():
            if not isinstance(key, str) or "." not in key:
                errors.append(
                    f"  {study_path}: parameter key '{key}' must use '<target>.<yaml.path>' format."
                )
                continue
            root = key.split(".", 1)[0]
            if root not in allowed_roots:
                errors.append(
                    f"  {study_path}: parameter key '{key}' must start with one of {sorted(allowed_roots)}."
                )
            if not isinstance(values, list) or len(values) == 0:
                errors.append(f"  {study_path}: parameters.{key} must be a non-empty list.")
 
    metrics = study_cfg.get("metrics", [])
    if metrics is not None and not isinstance(metrics, list):
        errors.append(f"  {study_path}: 'metrics' must be a list when provided.")
    elif isinstance(metrics, list):
        for i, metric in enumerate(metrics):
            if isinstance(metric, str):
                continue
            if not isinstance(metric, dict):
                errors.append(
                    f"  {study_path}: metrics[{i}] must be a string or mapping."
                )
                continue
            if "name" not in metric:
                errors.append(f"  {study_path}: metrics[{i}] missing required key 'name'.")
            if "source" not in metric:
                errors.append(f"  {study_path}: metrics[{i}] missing required key 'source'.")
 
    plotting = study_cfg.get("plotting", {})
    if plotting is not None and not isinstance(plotting, dict):
        errors.append(f"  {study_path}: 'plotting' must be a mapping when provided.")
    elif isinstance(plotting, dict):
        enabled = plotting.get("enabled")
        if enabled is not None and not isinstance(enabled, bool):
            errors.append(f"  {study_path}: plotting.enabled must be boolean when provided.")
        output_format = plotting.get("output_format")
        if output_format is not None and output_format not in {"png", "pdf", "svg"}:
            errors.append(f"  {study_path}: plotting.output_format must be one of ['png','pdf','svg'].")
 
    execution = study_cfg.get("execution", {})
    if execution is not None and not isinstance(execution, dict):
        errors.append(f"  {study_path}: 'execution' must be a mapping when provided.")
    elif isinstance(execution, dict):
        max_conc = execution.get("max_concurrent_array_tasks")
        if max_conc is not None and (not isinstance(max_conc, int) or max_conc <= 0):
            errors.append(
                f"  {study_path}: execution.max_concurrent_array_tasks must be a positive integer when provided."
            )
 
    if errors:
        _print_validation_errors(errors)
 
def _deep_set(container: dict, dotted_path: str, value):
    """Set nested dictionary value, creating intermediate maps when needed."""
    keys = dotted_path.split(".")
    current = container
    for key in keys[:-1]:
        if key not in current or not isinstance(current[key], dict):
            current[key] = {}
        current = current[key]
    current[keys[-1]] = value
 
def expand_parameter_matrix(parameters: dict) -> list:
    """Expand study parameter lists into cartesian-product combinations."""
    param_keys = list(parameters.keys())
    all_values = [parameters[k] for k in param_keys]
    combos = []
    for combo in itertools.product(*all_values):
        combos.append(dict(zip(param_keys, combo)))
    return combos
 
def get_cluster_total_tasks(cluster_cfg: dict) -> int:
    resources = cluster_cfg.get("resources", {})
    return int(resources.get("nodes", 1)) * int(resources.get("ntasks_per_node", 1))
 
def normalize_extension(ext: str) -> str:
    if ext is None:
        return None
    return str(ext).strip().lstrip(".")
 
def render_slurm_script(
    script_path: str,
    job_name: str,
    cluster_cfg: dict,
    command: list,
    workdir: str,
    stdout_path: str,
    stderr_path: str = None,
    env_vars: dict = None,
    array_spec: str = None
):
    """Render a Slurm batch script for a single command."""
    resources = cluster_cfg.get("resources", {})
    notifications = cluster_cfg.get("notifications", {}) or {}
    execution = cluster_cfg.get("execution", {}) or {}
    extra_sbatch = execution.get("extra_sbatch")
    module_setup = execution.get("module_setup", []) or []
 
    if stderr_path is None:
        stderr_path = stdout_path.replace(".out", ".err")
 
    lines = [
        "#!/bin/bash",
        f"#SBATCH --job-name={job_name}",
        f"#SBATCH --nodes={resources['nodes']}",
        f"#SBATCH --ntasks-per-node={resources['ntasks_per_node']}",
        f"#SBATCH --mem={resources['mem']}",
        f"#SBATCH --time={resources['time']}",
        f"#SBATCH --output={stdout_path}",
        f"#SBATCH --error={stderr_path}",
        f"#SBATCH --account={resources['account']}",
    ]
    partition = resources.get("partition")
    if partition:
        lines.append(f"#SBATCH --partition={partition}")
    if array_spec:
        lines.append(f"#SBATCH --array={array_spec}")
    mail_user = notifications.get("mail_user")
    mail_type = notifications.get("mail_type")
    if mail_user:
        lines.append(f"#SBATCH --mail-user={mail_user}")
    if mail_type:
        lines.append(f"#SBATCH --mail-type={mail_type}")
 
    if isinstance(extra_sbatch, dict):
        for key, value in extra_sbatch.items():
            flag = str(key)
            if not flag.startswith("--"):
                flag = f"--{flag}"
            if isinstance(value, bool):
                if value:
                    lines.append(f"#SBATCH {flag}")
            elif value is not None:
                lines.append(f"#SBATCH {flag}={value}")
    elif isinstance(extra_sbatch, list):
        for token in extra_sbatch:
            lines.append(f"#SBATCH {token}")
 
    lines.extend(
        [
            "",
            "set -euo pipefail",
            "",
            f"cd {shlex.quote(workdir)}",
            'echo "[$(date)] Starting job ${SLURM_JOB_NAME} (${SLURM_JOB_ID})"',
            'echo "[$(date)] Working directory: $PWD"',
        ]
    )
 
    for setup_line in module_setup:
        lines.append(str(setup_line))
 
    if env_vars:
        for key, value in env_vars.items():
            lines.append(f"export {key}={shlex.quote(str(value))}")
 
    cmd = " ".join(shlex.quote(str(tok)) for tok in command)
    lines.append(cmd)
    lines.append('echo "[$(date)] Job completed."')
 
    os.makedirs(os.path.dirname(script_path), exist_ok=True)
    with open(script_path, "w") as f:
        f.write("\n".join(lines) + "\n")
    os.chmod(script_path, 0o755)
 
def build_cluster_launch_command(cluster_cfg: dict, executable: str, executable_args: list) -> list:
    """Build scheduler launcher command (srun/mpirun/custom) from cluster config."""
    execution = cluster_cfg.get("execution", {}) or {}
    launcher = execution.get("launcher", "srun")
    launcher_args = [str(x) for x in execution.get("launcher_args", [])]
    ntasks = get_cluster_total_tasks(cluster_cfg)
 
    if launcher and launcher.lower() == "srun":
        has_n = any(token in {"-n", "--ntasks"} for token in launcher_args)
        cmd = ["srun"] + launcher_args
        if not has_n:
            cmd += ["-n", str(ntasks)]
        return cmd + [executable] + executable_args
 
    if launcher and launcher.lower() == "mpirun":
        has_np = any(token in {"-np", "-n"} for token in launcher_args)
        cmd = ["mpirun"] + launcher_args
        if not has_np:
            cmd += ["-np", str(ntasks)]
        return cmd + [executable] + executable_args
 
    # Custom launcher or no launcher.
    cmd = []
    if launcher:
        cmd.append(str(launcher))
    cmd += launcher_args
    cmd += [executable] + executable_args
    return cmd
 
def parse_slurm_job_id(sbatch_output: str) -> str:
    """Extract numeric job id from standard sbatch output."""
    match = re.search(r"Submitted batch job\s+(\d+)", sbatch_output or "")
    return match.group(1) if match else None
 
def submit_sbatch(script_path: str, dependency: str = None) -> dict:
    """Submit sbatch script and return submission metadata."""
    cmd = ["sbatch"]
    if dependency:
        cmd.append(f"--dependency=afterok:{dependency}")
    cmd.append(script_path)
    result = subprocess.run(cmd, text=True, capture_output=True, check=False)
    metadata = {
        "command": cmd,
        "returncode": result.returncode,
        "stdout": (result.stdout or "").strip(),
        "stderr": (result.stderr or "").strip(),
        "script": script_path,
    }
    if result.returncode != 0:
        print(f"[FATAL] sbatch submission failed for {script_path}\n{metadata['stderr']}", file=sys.stderr)
        sys.exit(result.returncode)
    metadata["job_id"] = parse_slurm_job_id(metadata["stdout"])
    if not metadata["job_id"]:
        print(
            f"[FATAL] Could not parse Slurm job id from sbatch output: {metadata['stdout']}",
            file=sys.stderr
        )
        sys.exit(1)
    return metadata
 
 
def _print_validation_errors(errors: list):
    """!
    @brief Prints validation errors and exits.
    @param[in] errors List of error message strings.
    """
    print(f"\n[FATAL] Configuration validation failed with {len(errors)} issue(s):", file=sys.stderr)
    for raw_error in errors:
        file_path, message = _split_error_file_and_message(raw_error)
        key_path = _extract_key_path(message)
        code = _classify_error_code(message)
        emit_structured_error(code, key=key_path, file_path=file_path, message=message)
    print(
        "\nHint: See examples/master_template/ for valid config structure and "
        "docs/pages/14_Config_Contract.md for key-level contract details.",
        file=sys.stderr,
    )
    sys.exit(1)
 
 
def generate_header(run_id: str, source_files: dict) -> str:
    """!
    @brief Creates a standard header block for all generated files.
    @param[in] run_id The unique identifier for the current simulation run.
    @param[in] source_files A dictionary of source profile files used.
    @return A formatted string containing the header.
    """
    header_parts = [
        "# ==============================================================================",
        "#                AUTO-GENERATED CONFIGURATION FILE",
        "# ------------------------------------------------------------------------------",
        f"#   Run ID:       {run_id}",
        f"#   Timestamp:    {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}",
        "#",
        "#   Source Configuration:"
    ]
    for name, path in source_files.items():
        header_parts.append(f"#     - {name:<12}: {os.path.basename(path)}")
    header_parts.extend([
        "#",
        "#   DO NOT EDIT THIS FILE MANUALLY. IT IS A MACHINE-READABLE ARTIFACT.",
        "# ==============================================================================\n"
    ])
    return "\n".join(header_parts)
 
def generate_simple_list_file(run_dir: str, run_id: str, cfg: dict, section: str, key: str, filename: str, header_sources: dict) -> str:
    """!
    @brief Generic function to create a file containing a simple list of strings.
    @param[in] run_dir The path to the main run directory.
    @param[in] run_id The unique identifier for the run.
    @param[in] cfg The dictionary containing the configuration data.
    @param[in] section The top-level key in the cfg dictionary.
    @param[in] key The second-level key whose value is the list of strings.
    @param[in] filename The name of the file to generate (e.g., 'whitelist.run').
    @param[in] header_sources A dictionary of source files for the header.
    @return The absolute path to the generated file.
    """
    print(f"[INFO] Generating {filename}...")
    config_dir = os.path.join(run_dir, "config")
    file_path = os.path.join(config_dir, filename)
    
    lines = [generate_header(run_id, header_sources)]
    items = cfg.get(section, {}).get(key, [])
    lines.extend(items)
 
    with open(file_path, "w") as f: f.write("\n".join(lines))
    print(f"[SUCCESS] Generated {filename}: {os.path.relpath(file_path)}")
    return os.path.abspath(file_path)
 
def generate_multi_block_bcs(run_dir: str, run_id: str, case_cfg: dict, source_files: dict) -> list:
    """!
    @brief Parses multi-block BCs from YAML, generates a .run file for each block,
           and returns a list of their absolute paths.
    @details Handles both simple list format (for single-block cases) and a
             list-of-lists (for multi-block cases) for boundary conditions.
    @param[in] run_dir The path to the main run directory.
    @param[in] run_id The unique identifier for the run.
    @param[in] case_cfg The parsed case.yml configuration dictionary.
    @param[in] source_files A dictionary of source files for the header.
    @return A list of absolute paths to the generated BC files.
    @throws ValueError if the number of BC definitions does not match the number of blocks.
    """
    print("[INFO] Generating boundary condition files...")
    config_dir = os.path.join(run_dir, "config")
    num_blocks = int(case_cfg.get('models', {}).get('domain', {}).get('blocks', 1))
    prepared_blocks = validate_and_prepare_boundary_conditions(case_cfg)
 
    generated_files = []
    for i, block_bcs_list in enumerate(prepared_blocks):
        file_name = "bcs.run" if num_blocks == 1 else f"bcs_block{i}.run"
        bcs_file_path = os.path.join(config_dir, file_name)
        bcs_lines = [generate_header(run_id, source_files)]
 
        for bc in block_bcs_list:
            face, bc_type, handler = bc['face'], bc['type'], bc['handler']
            params_str = ""
            if bc.get('params'):
                parts = []
                for k, v in bc['params'].items():
                    if isinstance(v, bool):
                        value_str = "true" if v else "false"
                    else:
                        value_str = str(v)
                    parts.append(f"{k}={value_str}")
                params_str = " ".join(parts)
            bcs_lines.append(f"{face:<20s} {bc_type:<12s} {handler:<20s} {params_str}")
        
        with open(bcs_file_path, "w") as f: f.write("\n".join(bcs_lines))
        
        print(f"[SUCCESS] Generated BCs for Block {i}: {os.path.relpath(bcs_file_path)}")
        generated_files.append(os.path.abspath(bcs_file_path))
        
    return generated_files
 
def format_flag_value(value):
    """!
    @brief Converts Python types to C-style command-line flag values.
    @param[in] value The Python object to convert (bool, list, or other).
    @return A string representation suitable for a C command-line parser.
    """
    if isinstance(value, bool):
        return "1" if value else "0"
    if isinstance(value, list):
        return ",".join(map(str, value))
    return str(value)
 
def normalize_momentum_solver_type(value: str) -> str:
    """!
    @brief Normalizes user-facing momentum solver names to C-enum CLI values.
    @param[in] value Human-readable or enum-like momentum solver string.
    @return Canonical value accepted by -mom_solver_type.
    @throws ValueError if the input cannot be mapped.
    """
    if value is None:
        raise ValueError("momentum solver type cannot be None")
 
    raw = str(value).strip()
    direct = {
        "EXPLICIT_RK",
        "DUALTIME_PICARD_RK4",
        "DUALTIME_NK_ARNOLDI",
        "DUALTIME_NK_ANALYTICAL_JACOBIAN"
    }
    if raw in direct:
        return raw
 
    normalized = raw.lower().replace("-", " ").replace("_", " ")
    normalized = " ".join(normalized.split())
    aliases = {
        "explicit rk4": "EXPLICIT_RK",
        "explicit rk": "EXPLICIT_RK",
        "explicit runge kutta 4": "EXPLICIT_RK",
        "explicit runge kutta": "EXPLICIT_RK",
        "dual time picard rk4": "DUALTIME_PICARD_RK4",
        "dualtime picard rk4": "DUALTIME_PICARD_RK4",
        "dual time picard": "DUALTIME_PICARD_RK4",
        "dualtime picard": "DUALTIME_PICARD_RK4",
        "dual time nk arnoldi": "DUALTIME_NK_ARNOLDI",
        "dualtime nk arnoldi": "DUALTIME_NK_ARNOLDI",
        "dual time nk analytical jacobian": "DUALTIME_NK_ANALYTICAL_JACOBIAN",
        "dualtime nk analytical jacobian": "DUALTIME_NK_ANALYTICAL_JACOBIAN"
    }
    mapped = aliases.get(normalized)
    if mapped is None:
        raise ValueError(
            f"Unknown momentum solver '{value}'. Use one of: "
            "'Explicit RK4', 'Dual Time Picard RK4', 'Dual Time NK Arnoldi', "
            "'Dual Time NK Analytical Jacobian' (or canonical enum values)."
        )
    return mapped
 
def normalize_field_init_mode(value: str) -> int:
    """!
    @brief Normalizes user-facing field init mode names to C enum/int codes (-finit).
    @param[in] value Human-readable or enum-like field initialization mode.
    @return Canonical integer code accepted by -finit.
    @throws ValueError if the input cannot be mapped.
    """
    if value is None:
        raise ValueError("field initialization mode cannot be None")
 
    raw = str(value).strip()
    direct = {"0": 0, "1": 1, "2": 2}
    if raw in direct:
        return direct[raw]
 
    normalized = raw.lower().replace("-", " ").replace("_", " ")
    normalized = " ".join(normalized.split())
    aliases = {
        "zero": 0,
        "constant": 1,
        "poiseuille": 2,
    }
    mapped = aliases.get(normalized)
    if mapped is None:
        raise ValueError(
            f"Unknown initial_conditions mode '{value}'. Use one of: 'Zero', 'Constant', 'Poiseuille' (or 0/1/2)."
        )
    return mapped
 
def normalize_particle_init_mode(value: str) -> int:
    """!
    @brief Normalizes particle init mode names to C enum/int codes (-pinit).
    @param[in] value Human-readable or enum-like particle initialization mode.
    @return Canonical integer code accepted by -pinit.
    @throws ValueError if the input cannot be mapped.
    """
    if value is None:
        raise ValueError("particle init mode cannot be None")
 
    raw = str(value).strip()
    direct = {"0": 0, "1": 1, "2": 2, "3": 3}
    if raw in direct:
        return direct[raw]
 
    normalized = raw.lower().replace("-", " ").replace("_", " ")
    normalized = " ".join(normalized.split())
    aliases = {
        "surface": 0,
        "surface random": 0,
        "volume": 1,
        "volumetric": 1,
        "point source": 2,
        "pointsource": 2,
        "surface edges": 3,
        "surfaceedges": 3,
        "surface at edges": 3,
    }
    mapped = aliases.get(normalized)
    if mapped is None:
        raise ValueError(
            f"Unknown particle init_mode '{value}'. Use one of: "
            "'Surface', 'Volume', 'PointSource', 'SurfaceEdges' (or 0/1/2/3)."
        )
    return mapped
 
def append_passthrough_flags(control_lines: list, options: dict):
    """!
    @brief Appends raw CLI flags to the control list from a {flag: value} dict.
    @details Boolean `true` is emitted as a switch with no value. Boolean `false`
             is skipped. All other values are emitted as "<flag> <value>".
    @param[out] control_lines The destination list of control-file lines.
    @param[in] options Mapping of raw CLI flags to values.
    """
    if not options:
        return
    for flag, value in options.items():
        if isinstance(value, bool):
            if value:
                control_lines.append(str(flag))
            continue
        control_lines.append(f"{flag} {format_flag_value(value)}")
 
def parse_and_add_model_flags(case_cfg: dict, control_lines: list):
    """!
    @brief Parses the 'models' section of case.yml and adds corresponding C-solver flags.
    @param[in] case_cfg The parsed case.yml configuration dictionary.
    @param[out] control_lines A list of strings to which C-flags will be appended.
    """
    models = case_cfg.get('models', {})
    FLAG_MAP = {
        'domain': {'blocks': '-nblk', 'i_periodic': '-i_periodic', 'j_periodic': '-j_periodic', 'k_periodic': '-k_periodic'},
        'physics.fsi': {'immersed': '-imm', 'moving_fsi': '-fsi'},
        'physics.particles': {'count': '-numParticles'},
        'physics.turbulence': {'les': '-les', 'rans': '-rans', 'wall_function': '-wallfunction'},
        'statistics': {'time_averaging': '-averaging'}
    }
    for section_path, flags in FLAG_MAP.items():
        current_level = models
        try:
            for key in section_path.split('.'): current_level = current_level[key]
            for yaml_key, flag in flags.items():
                if yaml_key in current_level:
                    control_lines.append(f"{flag} {format_flag_value(current_level[yaml_key])}")
        except KeyError: continue
 
    if models.get('physics', {}).get('dimensionality') == '2D':
        control_lines.append("-TwoD 1")
    
    particles_cfg = models.get('physics', {}).get('particles', {})
    p_init_mode_str = particles_cfg.get('init_mode', 'Surface')
    pinit_code = normalize_particle_init_mode(p_init_mode_str)
    control_lines.append(f"-pinit {pinit_code}")
    print(f"  - Particle Initialization Mode: {p_init_mode_str} (Code: {pinit_code})")
 
    if pinit_code == 2:
        point_cfg = particles_cfg.get('point_source', {})
        if not isinstance(point_cfg, dict):
            raise ValueError("models.physics.particles.point_source must be a mapping when init_mode is PointSource.")
        try:
            psrc_x = float(point_cfg['x'])
            psrc_y = float(point_cfg['y'])
            psrc_z = float(point_cfg['z'])
        except (KeyError, TypeError, ValueError):
            raise ValueError("PointSource init_mode requires numeric point_source.{x,y,z} values.")
        control_lines.append(f"-psrc_x {psrc_x}")
        control_lines.append(f"-psrc_y {psrc_y}")
        control_lines.append(f"-psrc_z {psrc_z}")
        print(f"  - Particle Point Source: ({psrc_x}, {psrc_y}, {psrc_z})")
 
    p_restart_mode = particles_cfg.get('restart_mode')
    if p_restart_mode:
        p_restart_mode_normalized = str(p_restart_mode).lower()
        if p_restart_mode_normalized not in {"init", "load"}:
            raise ValueError(f"Unknown particle restart_mode '{p_restart_mode}'. Options are 'init' or 'load'.")
        control_lines.append(f"-particle_restart_mode \"{p_restart_mode}\"")
 
def parse_solver_config(solver_cfg: dict) -> dict:
    """!
    @brief Parses the structured solver.yml into a flat dictionary of {flag: value}.
    @param[in] solver_cfg The parsed solver.yml configuration dictionary.
    @return A dictionary where keys are C-solver flags and values are the corresponding settings.
    """
    flags = {}
    if 'operation_mode' in solver_cfg and isinstance(solver_cfg['operation_mode'], dict):
        op_mode = solver_cfg['operation_mode']
        if 'eulerian_field_source' in op_mode:
            flags['-euler_field_source'] = f"\"{op_mode.get('eulerian_field_source', 'solve')}\""
        if 'analytical_type' in op_mode and op_mode.get('analytical_type') is not None:
            flags['-analytical_type'] = f"\"{op_mode.get('analytical_type')}\""
 
    selected_solver = None
    if 'strategy' in solver_cfg:
        s = solver_cfg['strategy']
        if 'central_diff' in s:
            flags['-central'] = format_flag_value(s['central_diff'])
        # Preferred selector.
        if 'momentum_solver' in s:
            selected_solver = normalize_momentum_solver_type(s['momentum_solver'])
        elif 'implicit' in s:
            raise ValueError("Legacy key 'strategy.implicit' is not supported. Use 'strategy.momentum_solver'.")
 
    ms = solver_cfg.get('momentum_solver', {})
    if isinstance(ms, dict) and 'type' in ms and selected_solver is None:
        selected_solver = normalize_momentum_solver_type(ms['type'])
 
    if selected_solver is None:
        selected_solver = "DUALTIME_PICARD_RK4"
    flags['-mom_solver_type'] = f"\"{selected_solver}\""
 
    if 'tolerances' in solver_cfg:
        t = solver_cfg['tolerances']
        tol_map = {
            'max_iterations': '-mom_max_pseudo_steps',
            'absolute_tol': '-mom_atol',
            'relative_tol': '-mom_rtol',
            'step_tol': '-imp_stol'
        }
        for key, flag in tol_map.items():
            if key in t:
                flags[flag] = t[key]
 
    def _append_dualtime_options(cfg: dict):
        if 'max_pseudo_steps' in cfg:
            flags['-mom_max_pseudo_steps'] = cfg['max_pseudo_steps']
        if 'absolute_tol' in cfg:
            flags['-mom_atol'] = cfg['absolute_tol']
        if 'relative_tol' in cfg:
            flags['-mom_rtol'] = cfg['relative_tol']
        if 'step_tol' in cfg:
            flags['-imp_stol'] = cfg['step_tol']
        if 'pseudo_cfl' in cfg:
            pcfl = cfg['pseudo_cfl']
            if 'initial' in pcfl:
                flags['-pseudo_cfl'] = pcfl['initial']
            if 'minimum' in pcfl:
                flags['-min_pseudo_cfl'] = pcfl['minimum']
            if 'maximum' in pcfl:
                flags['-max_pseudo_cfl'] = pcfl['maximum']
            if 'growth_factor' in pcfl:
                flags['-pseudo_cfl_growth_factor'] = pcfl['growth_factor']
            if 'reduction_factor' in pcfl:
                flags['-pseudo_cfl_reduction_factor'] = pcfl['reduction_factor']
        if 'rk4_residual_noise_allowance_factor' in cfg:
            flags['-mom_dt_rk4_residual_norm_noise_allowance_factor'] = cfg['rk4_residual_noise_allowance_factor']
 
    if isinstance(ms, dict):
        allowed_ms_keys = {'type', 'dual_time_picard_rk4'}
        unknown_ms_keys = sorted(set(ms.keys()) - allowed_ms_keys)
        if unknown_ms_keys:
            raise ValueError(
                f"Unsupported momentum_solver keys/blocks: {unknown_ms_keys}. "
                "Currently supported block: 'dual_time_picard_rk4'."
            )
 
        dt_picard_cfg = ms.get('dual_time_picard_rk4')
        if dt_picard_cfg is not None:
            if selected_solver != "DUALTIME_PICARD_RK4":
                raise ValueError(
                    f"momentum_solver.dual_time_picard_rk4 is set but selected solver is {selected_solver}."
                )
            if not isinstance(dt_picard_cfg, dict):
                raise ValueError("momentum_solver.dual_time_picard_rk4 must be a mapping.")
            _append_dualtime_options(dt_picard_cfg)
    if 'pressure_solver' in solver_cfg:
        ps = solver_cfg['pressure_solver']
        if 'tolerance' in ps: flags['-poisson_tol'] = ps['tolerance']
        if 'multigrid' in ps:
            mg = ps['multigrid']
            mg_map = {'levels': '-mg_level', 'pre_sweeps': '-mg_pre_it', 'post_sweeps': '-mg_post_it'}
            for key, flag in mg_map.items():
                if key in mg: flags[flag] = mg[key]
            if 'semi_coarsening' in mg:
                sc = mg['semi_coarsening']
                if 'i' in sc: flags['-mg_i_semi'] = format_flag_value(sc['i'])
                if 'j' in sc: flags['-mg_j_semi'] = format_flag_value(sc['j'])
                if 'k' in sc: flags['-mg_k_semi'] = format_flag_value(sc['k'])
            if 'level_solvers' in mg:
                for level_name, settings in mg['level_solvers'].items():
                    level_num = level_name.split('_')[-1]
                    for key, value in settings.items():
                        flags[f"-ps_mg_levels_{level_num}_{key}"] = format_flag_value(value)
    if 'petsc_passthrough_options' in solver_cfg:
        passthrough = solver_cfg['petsc_passthrough_options']
        if passthrough:
            for key, value in passthrough.items():
                flags[key] = format_flag_value(value)
    return flags
 
def generate_solver_control_file(run_dir, run_id, configs, num_procs, monitor_files):
    """!
    @brief Generates the main .control file for the C-solver.
    @details Orchestrates the conversion of all YAML configurations (case, solver, monitor)
             into a single, machine-readable file of command-line flags.
    @param[in] run_dir The path to the main run directory.
    @param[in] run_id The unique identifier for the run.
    @param[in] configs A dictionary containing the parsed YAML data.
    @param[in] num_procs The number of MPI processes for the run.
    @param[in] monitor_files A dictionary containing paths to generated monitor files.
    @return The absolute path to the generated solver control file.
    """
    print("[INFO] Generating master solver control file...")
    case_cfg, solver_cfg, monitor_cfg = configs['case'], configs['solver'], configs['monitor']
    source_files = {'Case': configs['case_path'], 'Solver': configs['solver_path'], 'Monitor': configs['monitor_path']}
    
    control_lines = []
    try:
        props, run_ctrl = case_cfg['properties'], case_cfg['run_control']
        scales, fluid, ic = props['scaling'], props['fluid'], props['initial_conditions']
        L_ref, U_ref, rho, mu = float(scales['length_ref']), float(scales['velocity_ref']), float(fluid['density']), float(fluid['viscosity'])
        reynolds = (rho * U_ref * L_ref) / mu if mu != 0 else float('inf')
        dt_phys = float(run_ctrl['dt_physical'])
        T_ref = L_ref / U_ref if U_ref != 0 else float('inf')
        dt_nondim = dt_phys / T_ref if T_ref != float('inf') else 0.0
        u, v, w = float(ic['u_physical']), float(ic['v_physical']), float(ic['w_physical'])
        finit_mode_str = ic.get('mode', 'Constant')
        finit_code = normalize_field_init_mode(finit_mode_str)
        print(f"  - Reynolds Number (Re) = {reynolds:.4f}")
        print(f"  - Non-Dimensional dt*  = {dt_nondim:.6f}")
        print(f"  - Field Initialization Mode: {finit_mode_str} (Code: {finit_code})")
        control_lines.extend([
            f"-start_step {run_ctrl['start_step']}", f"-totalsteps {run_ctrl['total_steps']}",
            f"-ren {reynolds}", f"-dt {dt_nondim}", f"-finit {finit_code}",
            f"-ucont_x {u/U_ref if U_ref!=0 else 0}", f"-ucont_y {v/U_ref if U_ref!=0 else 0}", f"-ucont_z {w/U_ref if U_ref!=0 else 0}",
            f"-scaling_L_ref {L_ref}", f"-scaling_U_ref {U_ref}", f"-scaling_rho_ref {rho}"
        ])
    except (KeyError, TypeError, ZeroDivisionError, ValueError) as e:
        print(f"[FATAL] Error processing case.yml: {e}", file=sys.stderr)
        sys.exit(1)
        
    try:
        bcs_files = generate_multi_block_bcs(run_dir, run_id, case_cfg, source_files)
    except ValueError as e:
        print(f"[FATAL] Invalid boundary_conditions in case.yml: {e}", file=sys.stderr)
        sys.exit(1)
    control_lines.append(f"-bcs_files \"{','.join(bcs_files)}\"")
 
    # --- CORRECTED: Add paths for whitelist and profile files ---
    control_lines.append(f"-whitelist_config_file {monitor_files['whitelist']}")
    control_lines.append(f"-profile_config_file {monitor_files['profile']}")
    
    grid_cfg = case_cfg.get('grid', {})
    grid_mode = grid_cfg.get('mode')
    expected_nblk = int(case_cfg.get('models', {}).get('domain', {}).get('blocks', 1))
 
    if grid_mode == 'file':
        print("[INFO] Grid Mode: Using external file...")
        case_file_dir = os.path.dirname(configs['case_path'])
        source_grid = grid_cfg['source_file']
        if not os.path.isabs(source_grid):
            source_grid = os.path.abspath(os.path.join(case_file_dir, source_grid))
        nondim_grid_path = os.path.join(run_dir, "config", "grid.run")
        try:
            summary = validate_and_nondimensionalize_picgrid(
                source_grid, nondim_grid_path, L_ref, expected_nblk=expected_nblk
            )
            print(
                f"[SUCCESS] Validated and non-dimensionalized grid: {os.path.relpath(nondim_grid_path)} "
                f"(nblk={summary['nblk']}, total_nodes={summary['total_nodes']})"
            )
            control_lines.append(f"-grid_file {nondim_grid_path}")
        except Exception as e:
            print(f"[FATAL] Failed to process grid file '{source_grid}': {e}", file=sys.stderr)
            sys.exit(1)
    elif grid_mode == 'grid_gen':
        print("[INFO] Grid Mode: Generating external grid via grid.gen...")
        nondim_grid_path = os.path.join(run_dir, "config", "grid.run")
        try:
            generated_grid = run_grid_generator(configs['case_path'], run_dir, grid_cfg)
            summary = validate_and_nondimensionalize_picgrid(
                generated_grid, nondim_grid_path, L_ref, expected_nblk=expected_nblk
            )
            print(
                f"[SUCCESS] grid.gen output validated and non-dimensionalized: {os.path.relpath(nondim_grid_path)} "
                f"(nblk={summary['nblk']}, total_nodes={summary['total_nodes']})"
            )
            control_lines.append(f"-grid_file {nondim_grid_path}")
        except Exception as e:
            print(f"[FATAL] Grid generation failed: {e}", file=sys.stderr)
            sys.exit(1)
    elif grid_mode == 'programmatic_c':
        print("[INFO] Grid Mode: Programmatic C...")
        grid_settings = dict(grid_cfg.get('programmatic_settings', {}))
        control_lines.append("-grid")
        px, py, pz = grid_settings.get('da_processors_x'), grid_settings.get('da_processors_y'), grid_settings.get('da_processors_z')
        if any(isinstance(p, (list, tuple)) for p in [px, py, pz] if p is not None):
            raise ValueError(
                "da_processors_x/y/z must be scalar integers. "
                "Per-block MPI decomposition is not implemented on the C side; DMDA layout is global."
            )
        if num_procs > 1 and all(p is not None for p in [px, py, pz]):
            if not all(isinstance(p, int) and p > 0 for p in [px, py, pz]):
                raise ValueError("da_processors_x/y/z must be positive integers when provided.")
            total_layout = px * py * pz
            if total_layout != num_procs:
                raise ValueError(f"Processor layout mismatch: product ({total_layout}) != processes ({num_procs}).")
            print(f"[INFO] Applying user-defined processor layout for {num_procs} processes.")
        else:
            if num_procs == 1: print("[INFO] Serial run, ignoring da_processors layout.")
            else: print("[INFO] Letting PETSc automatically determine processor layout.")
            for p_key in ['da_processors_x', 'da_processors_y', 'da_processors_z']:
                grid_settings.pop(p_key, None)
        for key, value in grid_settings.items(): control_lines.append(f"-{key} {format_flag_value(value)}")
    else:
        raise ValueError(f"Unknown or missing grid mode '{grid_mode}' in case.yml.")
    
    parse_and_add_model_flags(case_cfg, control_lines)
    
    if 'solver_parameters' in case_cfg:
        params = case_cfg['solver_parameters']
        if params:
            for key, value in params.items():
                control_lines.append(f"{key} {format_flag_value(value)}")
 
    try:
        solver_flags = parse_solver_config(solver_cfg)
    except ValueError as e:
        print(f"[FATAL] Invalid solver.yml settings: {e}", file=sys.stderr)
        sys.exit(1)
    for flag, value in solver_flags.items(): control_lines.append(f"{flag} {value}")
 
    append_passthrough_flags(control_lines, monitor_cfg.get('solver_monitoring', {}))
    
    io_cfg = monitor_cfg.get('io', {})
    if 'data_output_frequency' in io_cfg: control_lines.append(f"-tio {io_cfg['data_output_frequency']}")
    if 'particle_log_interval' in io_cfg: control_lines.append(f"-logfreq {io_cfg['particle_log_interval']}")
    if 'directories' in io_cfg:
        dirs = io_cfg['directories']
        if 'output' in dirs: control_lines.append(f"-output_dir {dirs['output']}")
        if 'restart' in dirs: control_lines.append(f"-restart_dir {dirs['restart']}")
        if 'log' in dirs: control_lines.append(f"-log_dir {dirs['log']}")
        if 'eulerian_subdir' in dirs: control_lines.append(f"-euler_subdir {dirs['eulerian_subdir']}")
        if 'particle_subdir' in dirs: control_lines.append(f"-particle_subdir {dirs['particle_subdir']}")
        
    final_content = generate_header(run_id, source_files) + "\n".join(control_lines)
    control_file_path = os.path.join(run_dir, "config", f"{run_id}.control")
    with open(control_file_path, "w") as f: f.write(final_content)
    print(f"[SUCCESS] Generated solver control file: {os.path.relpath(control_file_path)}")
    return os.path.abspath(control_file_path)
 
def generate_post_recipe_file(run_dir: str, run_id: str, post_cfg: dict, source_files: dict) -> str:
    """!
    @brief Generates a key=value config file (post.run) for the C post-processor.
    @details Translates the structured post-processing YAML into the specific flat
             key-value format required by the C executable, including complex,
             semicolon-separated pipeline strings.
    @param[in] run_dir The path to the main run directory.
    @param[in] run_id The unique identifier for the run.
    @param[in] post_cfg The parsed post-profile YAML configuration dictionary.
    @param[in] source_files A dictionary of source files for the header.
    @return The absolute path to the generated post.run recipe file.
    """
    print("[INFO] Generating post-processor recipe file (post.run)...")
    config_dir = os.path.join(run_dir, "config")
    post_recipe_path = os.path.join(config_dir, "post.run")
    
    lines = [generate_header(run_id, source_files)]
 
    c_config = {}
 
    # --- 1. Process Run Control ---
    # Accept snake_case names (preferred) with camelCase fallback for backwards compatibility.
    rc = post_cfg.get('run_control', {})
    c_config['startTime'] = rc.get('start_step', rc.get('startTime', 0))
    c_config['endTime'] = rc.get('end_step', rc.get('endTime', 10))
    c_config['timeStep'] = rc.get('step_interval', rc.get('timeStep', 1))
 
    # --- 2. Process Global Operations ---
    eulerian_pipeline_parts = []
    if post_cfg.get('global_operations', {}).get('dimensionalize', False):
        eulerian_pipeline_parts.append('DimensionalizeAllLoadedFields')
 
    # --- 3. Build Eulerian Pipeline String ---
    for task in post_cfg.get('eulerian_pipeline', []):
        task_name = task.get('task')
        if task_name == 'q_criterion':
            eulerian_pipeline_parts.append('ComputeQCriterion')
        elif task_name == 'normalize_field':
            field = task.get('field', 'P')
            eulerian_pipeline_parts.append(f'NormalizeRelativeField:{field}')
            ref_point = task.get('reference_point', [1, 1, 1])
            c_config['reference_ip'] = ref_point[0]
            c_config['reference_jp'] = ref_point[1]
            c_config['reference_kp'] = ref_point[2]
        elif task_name == 'nodal_average':
            in_field = task.get('input_field')
            out_field = task.get('output_field')
            if in_field and out_field:
                eulerian_pipeline_parts.append(f'CellToNodeAverage:{in_field}>{out_field}')
 
    if eulerian_pipeline_parts:
        c_config['process_pipeline'] = ";".join(eulerian_pipeline_parts)
 
    # --- 4. Build Lagrangian Pipeline String ---
    lagrangian_pipeline_parts = []
    for task in post_cfg.get('lagrangian_pipeline', []):
        task_name = task.get('task')
        if task_name == 'specific_ke':
            in_field = task.get('input_field')
            out_field = task.get('output_field')
            if in_field and out_field:
                lagrangian_pipeline_parts.append(f'ComputeSpecificKE:{in_field}>{out_field}')
    
    if lagrangian_pipeline_parts:
        c_config['particle_pipeline'] = ";".join(lagrangian_pipeline_parts)
 
    # --- 4B. Build Statistics Pipeline String ---
    statistics_pipeline_parts = []
    statistics_output_prefix = None
    stats_cfg = post_cfg.get('statistics_pipeline')
    stats_entries = []
    if isinstance(stats_cfg, list):
        stats_entries = stats_cfg
    elif isinstance(stats_cfg, dict):
        stats_entries = stats_cfg.get('tasks', [])
        statistics_output_prefix = stats_cfg.get('output_prefix')
 
    for entry in stats_entries:
        if isinstance(entry, str):
            task_name = entry
        elif isinstance(entry, dict):
            task_name = entry.get('task')
        else:
            continue
        try:
            statistics_pipeline_parts.append(normalize_statistics_task(task_name))
        except ValueError:
            # validation should catch this earlier; keep generation tolerant.
            continue
 
    if statistics_pipeline_parts:
        c_config['statistics_pipeline'] = ";".join(statistics_pipeline_parts)
    if statistics_output_prefix is None:
        statistics_output_prefix = post_cfg.get('statistics_output_prefix')
    if statistics_output_prefix:
        c_config['statistics_output_prefix'] = statistics_output_prefix
 
    # --- 5. Process I/O ---
    io = post_cfg.get('io', {})
    c_config['output_prefix'] = io.get('output_directory','viz')+'/'+io.get('output_filename_prefix', 'Field')
    c_config['particle_output_prefix'] = io.get('output_directory','viz')+'/'+io.get('particle_filename_prefix', 'Particle')
    c_config['output_particles'] = io.get('output_particles', False)
    c_config['particle_output_freq'] = io.get('particle_subsampling_frequency', 1)
    c_config['output_fields_instantaneous'] = ",".join(io.get('eulerian_fields', []))
    c_config['output_fields_averaged'] = ",".join(io.get('eulerian_fields_averaged', []))
    c_config['particle_fields_instantaneous'] = ",".join(io.get('particle_fields', []))
    input_extensions = io.get('input_extensions', {})
    if isinstance(input_extensions, dict):
        e_ext = input_extensions.get('eulerian')
        p_ext = input_extensions.get('particle')
        if e_ext:
            c_config['eulerianExt'] = str(e_ext).strip().lstrip('.')
        if p_ext:
            c_config['particleExt'] = str(p_ext).strip().lstrip('.')
 
    # --- 6. Add Source Directory ---
    if 'source_data' in post_cfg and 'directory' in post_cfg['source_data']:
        c_config['source_directory'] = post_cfg['source_data']['directory']
 
    # --- 7. Write the final file ---
    for key, value in c_config.items():
        if value is not None and str(value) != "":
            lines.append(f"{key} = {value}")
 
    with open(post_recipe_path, "w") as f: f.write("\n".join(lines))
    print(f"[SUCCESS] Generated post-processor recipe: {os.path.relpath(post_recipe_path)}")
    return os.path.abspath(post_recipe_path)
 
def execute_command(command: list, run_dir: str, log_filename: str, monitor_cfg: dict = None):
    """!
    @brief Executes a command, streaming its output to the console and a log file.
    ...
    @param[in] monitor_cfg Optional. If provided, used to set LOG_LEVEL in a custom environment.
                           If None, the process inherits the parent's environment directly.
    """
    # Create the log directory if it doesn't exist.
    log_dir = os.path.join(run_dir, "logs")
    os.makedirs(log_dir, exist_ok=True)
    
    log_path = os.path.join(log_dir, log_filename)
    print(f"[INFO] Launching Command...\n  > {' '.join(command)}")
    print(f"       Log file: {os.path.relpath(log_path)}")
    print("-" * 60)
 
    # --- Environment Handling ---
    popen_kwargs = {
        "stdout": subprocess.PIPE, "stderr": subprocess.STDOUT,
        "cwd": run_dir, "bufsize": 1, "universal_newlines": True,
        "encoding": 'utf-8', "errors": 'replace'
    }
 
    if monitor_cfg:
        print("[INFO] Creating custom environment to set LOG_LEVEL.")
        run_env = os.environ.copy()
        verbosity = monitor_cfg.get('logging', {}).get('verbosity', 'INFO').upper()
        run_env['LOG_LEVEL'] = verbosity
        print(f"[INFO] Setting LOG_LEVEL={verbosity} for C executable.")
        popen_kwargs['env'] = run_env
    else:
        print("[INFO] Using inherited environment for process.")
 
    print("-" * 60)
    try:
        # Pass the constructed keyword arguments dictionary to Popen
        process = subprocess.Popen(command, **popen_kwargs)
        
        with open(log_path, "w") as log_file:
            for line in process.stdout:
                sys.stdout.write(line)
                log_file.write(line)
        process.wait()
        return_code = process.returncode
        print("-" * 60)
        if return_code == 0:
            print(f"[SUCCESS] Execution finished successfully.")
        else:
            print(f"[FATAL] Execution failed with exit code {return_code}. Check log: {os.path.relpath(log_path)}", file=sys.stderr)
            sys.exit(return_code)
    except FileNotFoundError:
        print(f"[FATAL] Command not found or is not executable: '{command[0]}'", file=sys.stderr)
        print("        Please check that the path is correct and the file has execute permissions.", file=sys.stderr)
        sys.exit(1)
    except Exception as e:
        print(f"[FATAL] An unexpected error occurred during execution: {e}", file=sys.stderr)
        sys.exit(1)
 
def auto_identify_run_inputs(config_dir: str):
    """Auto-detect case.yml, monitor.yml, and *.control in a run config directory."""
    all_yml_files = glob.glob(os.path.join(config_dir, "*.yml"))
    case_path, monitor_path = None, None
    for f_path in all_yml_files:
        try:
            content = read_yaml_file(f_path)
            if not isinstance(content, dict):
                continue
            if 'models' in content and 'boundary_conditions' in content:
                case_path = f_path
            elif 'io' in content and 'logging' in content:
                monitor_path = f_path
        except Exception as e:
            print(f"[WARNING] Could not parse or inspect '{f_path}': {e}", file=sys.stderr)
    try:
        solver_control_path = glob.glob(os.path.join(config_dir, "*.control"))[0]
    except IndexError:
        solver_control_path = None
    return case_path, monitor_path, solver_control_path
 
def resolve_post_source_directory(run_dir: str, monitor_cfg: dict, post_cfg: dict, strict: bool = True) -> str:
    """Resolve post source directory token and optionally enforce existence."""
    solver_output_dir_rel = monitor_cfg.get('io', {}).get('directories', {}).get('output', 'results')
    solver_output_dir_abs = os.path.join(run_dir, solver_output_dir_rel)
    source_dir_template = post_cfg.get('source_data', {}).get('directory', '<solver_output_dir>')
    if source_dir_template == '<solver_output_dir>':
        resolved_source_dir = solver_output_dir_abs
        print(f"[INFO] Post-processor source data: {os.path.relpath(resolved_source_dir)}")
    else:
        resolved_source_dir = os.path.abspath(os.path.join(run_dir, source_dir_template))
        print(f"[INFO] Post-processor source data (user-defined): {os.path.relpath(resolved_source_dir)}")
 
    if strict and (not os.path.isdir(resolved_source_dir) or not os.listdir(resolved_source_dir)):
        print(
            f"[FATAL] Source data directory for post-processing not found or empty: {os.path.relpath(resolved_source_dir)}",
            file=sys.stderr
        )
        sys.exit(1)
    if not strict and (not os.path.isdir(resolved_source_dir) or not os.listdir(resolved_source_dir)):
        print("[WARNING] Source data directory is not available yet; keeping deferred path for scheduled post job.")
    return resolved_source_dir
 
def render_slurm_array_stage_script(
    script_path: str,
    job_name: str,
    cluster_cfg: dict,
    array_spec: str,
    case_index_tsv: str,
    stage: str,
    solver_exe: str,
    post_exe: str,
    stdout_path: str,
    stderr_path: str
):
    """Render array script that maps SLURM_ARRAY_TASK_ID to per-case run artifacts."""
    resources = cluster_cfg.get("resources", {})
    notifications = cluster_cfg.get("notifications", {}) or {}
    execution = cluster_cfg.get("execution", {}) or {}
    module_setup = execution.get("module_setup", []) or []
    extra_sbatch = execution.get("extra_sbatch")
 
    lines = [
        "#!/bin/bash",
        f"#SBATCH --job-name={job_name}",
        f"#SBATCH --nodes={resources['nodes']}",
        f"#SBATCH --ntasks-per-node={resources['ntasks_per_node']}",
        f"#SBATCH --mem={resources['mem']}",
        f"#SBATCH --time={resources['time']}",
        f"#SBATCH --output={stdout_path}",
        f"#SBATCH --error={stderr_path}",
        f"#SBATCH --account={resources['account']}",
        f"#SBATCH --array={array_spec}",
    ]
    partition = resources.get("partition")
    if partition:
        lines.append(f"#SBATCH --partition={partition}")
    mail_user = notifications.get("mail_user")
    mail_type = notifications.get("mail_type")
    if mail_user:
        lines.append(f"#SBATCH --mail-user={mail_user}")
    if mail_type:
        lines.append(f"#SBATCH --mail-type={mail_type}")
    if isinstance(extra_sbatch, dict):
        for key, value in extra_sbatch.items():
            flag = str(key)
            if not flag.startswith("--"):
                flag = f"--{flag}"
            if isinstance(value, bool):
                if value:
                    lines.append(f"#SBATCH {flag}")
            elif value is not None:
                lines.append(f"#SBATCH {flag}={value}")
    elif isinstance(extra_sbatch, list):
        for token in extra_sbatch:
            lines.append(f"#SBATCH {token}")
 
    lines.extend([
        "",
        "set -euo pipefail",
        "",
        f'CASE_INDEX_FILE={shlex.quote(case_index_tsv)}',
        'LINE=$(sed -n "$((SLURM_ARRAY_TASK_ID + 1))p" "$CASE_INDEX_FILE")',
        'if [ -z "$LINE" ]; then',
        '  echo "No case entry for array index ${SLURM_ARRAY_TASK_ID}" >&2',
        '  exit 1',
        "fi",
        "IFS=$'\\t' read -r CASE_INDEX CASE_ID RUN_DIR CONTROL_FILE POST_RECIPE_FILE LOG_LEVEL POST_PREFIX <<< \"$LINE\"",
        'cd "$RUN_DIR"',
        'echo "[$(date)] Starting case ${CASE_ID} (array index ${SLURM_ARRAY_TASK_ID})"',
        'export LOG_LEVEL="${LOG_LEVEL}"',
    ])
 
    for setup_line in module_setup:
        lines.append(str(setup_line))
 
    if stage == "solve":
        cmd = build_cluster_launch_command(
            cluster_cfg,
            solver_exe,
            ["-control_file", "$CONTROL_FILE"]
        )
    else:
        cmd = build_cluster_launch_command(
            cluster_cfg,
            post_exe,
            ["-control_file", "$CONTROL_FILE", "-postprocessing_config_file", "$POST_RECIPE_FILE"]
        )
 
    # Keep shell variables unresolved inside sbatch script.
    def _token(tok: str) -> str:
        if tok.startswith("$"):
            return tok
        return shlex.quote(str(tok))
 
    lines.append(" ".join(_token(t) for t in cmd))
    lines.append('echo "[$(date)] Completed case ${CASE_ID}"')
 
    os.makedirs(os.path.dirname(script_path), exist_ok=True)
    with open(script_path, "w") as f:
        f.write("\n".join(lines) + "\n")
    os.chmod(script_path, 0o755)
 
def extract_metric_from_csv(case_dir: str, spec: dict):
    """Extract a scalar metric from a CSV source."""
    file_glob = spec.get("file_glob", "**/*_msd.csv")
    candidates = sorted(glob.glob(os.path.join(case_dir, file_glob), recursive=True))
    if not candidates:
        return None
    csv_path = candidates[0]
    rows = []
    with open(csv_path, "r", newline="") as f:
        reader = csv.DictReader(f)
        if reader.fieldnames:
            for row in reader:
                rows.append(row)
            if not rows:
                return None
            column = spec.get("column")
            if not column:
                for name in reversed(reader.fieldnames):
                    if name and name.lower() not in {"step", "time", "timestep"}:
                        column = name
                        break
            if not column:
                return None
            values = []
            for row in rows:
                try:
                    values.append(float(row[column]))
                except Exception:
                    continue
        else:
            return None
    if not values:
        return None
    reduction = str(spec.get("reduction", "last")).lower()
    if reduction == "mean":
        return float(np.mean(values))
    if reduction == "min":
        return float(np.min(values))
    if reduction == "max":
        return float(np.max(values))
    return float(values[-1])
 
def extract_metric_from_log(case_dir: str, spec: dict):
    """Extract a scalar metric from a log file using regex."""
    file_glob = spec.get("file_glob", "logs/*.log")
    regex = spec.get("regex")
    if not regex:
        return None
    candidates = sorted(glob.glob(os.path.join(case_dir, file_glob), recursive=True))
    if not candidates:
        return None
    pattern = re.compile(regex)
    values = []
    for path in candidates:
        try:
            with open(path, "r", encoding="utf-8", errors="replace") as f:
                for line in f:
                    m = pattern.search(line)
                    if m:
                        try:
                            values.append(float(m.group(1)))
                        except Exception:
                            pass
        except OSError:
            continue
    if not values:
        return None
    reduction = str(spec.get("reduction", "last")).lower()
    if reduction == "mean":
        return float(np.mean(values))
    if reduction == "min":
        return float(np.min(values))
    if reduction == "max":
        return float(np.max(values))
    return float(values[-1])
 
def normalize_metric_spec(metric):
    """Normalize study metric definitions to a common dictionary form."""
    if isinstance(metric, str):
        if metric.lower() in {"msd", "msd_final"}:
            return {
                "name": "msd_final",
                "source": "statistics_csv",
                "file_glob": "**/*_msd.csv",
                "reduction": "last",
            }
        return {"name": metric, "source": "log_regex", "regex": metric}
    return dict(metric)
 
def aggregate_study_metrics(study_cfg: dict, cases: list, results_dir: str) -> str:
    """Collect metric values from generated case directories into one CSV."""
    metrics = study_cfg.get("metrics", [])
    if not metrics:
        metrics = ["msd_final"]
    normalized_specs = [normalize_metric_spec(m) for m in metrics]
 
    rows = []
    for case in cases:
        row = {"case_id": case["case_id"]}
        for p_key, p_val in case["parameters"].items():
            row[p_key] = p_val
        for spec in normalized_specs:
            name = spec.get("name", "metric")
            source = str(spec.get("source", "")).lower()
            if source in {"statistics_csv", "csv"}:
                row[name] = extract_metric_from_csv(case["run_dir"], spec)
            elif source in {"log_regex", "log"}:
                row[name] = extract_metric_from_log(case["run_dir"], spec)
            else:
                row[name] = None
        rows.append(row)
 
    if not rows:
        return None
 
    all_keys = []
    seen = set()
    for row in rows:
        for k in row.keys():
            if k not in seen:
                seen.add(k)
                all_keys.append(k)
 
    os.makedirs(results_dir, exist_ok=True)
    out_csv = os.path.join(results_dir, "metrics_table.csv")
    with open(out_csv, "w", newline="") as f:
        writer = csv.DictWriter(f, fieldnames=all_keys)
        writer.writeheader()
        writer.writerows(rows)
    print(f"[SUCCESS] Aggregated metrics table: {os.path.relpath(out_csv)}")
    return out_csv
 
def infer_plot_x_axis(study_cfg: dict, rows: list):
    """Infer x-axis key/values for study plots."""
    params = list((study_cfg.get("parameters") or {}).keys())
    if not params or not rows:
        return None, None
 
    study_type = study_cfg.get("study_type")
    if study_type == "grid_independence":
        has_im = "case.grid.programmatic_settings.im" in params
        has_jm = "case.grid.programmatic_settings.jm" in params
        has_km = "case.grid.programmatic_settings.km" in params
        if has_im and has_jm and has_km:
            xs = []
            for row in rows:
                try:
                    im = float(row["case.grid.programmatic_settings.im"])
                    jm = float(row["case.grid.programmatic_settings.jm"])
                    km = float(row["case.grid.programmatic_settings.km"])
                    xs.append((im * jm * km) ** (1.0 / 3.0))
                except Exception:
                    return None, None
            return "N^(1/3)", xs
 
    primary = params[0]
    xs = []
    for row in rows:
        try:
            xs.append(float(row[primary]))
        except Exception:
            return None, None
    return primary, xs
 
def generate_study_plots(study_cfg: dict, metrics_csv: str, plots_dir: str):
    """Generate metric-vs-parameter plots for completed studies."""
    plotting_cfg = study_cfg.get("plotting", {}) or {}
    if plotting_cfg.get("enabled", True) is False:
        print("[INFO] Plotting disabled by study.yml.")
        return []
    if plt is None:
        print("[WARNING] matplotlib not available; skipping plot generation.")
        return []
    if not metrics_csv or not os.path.isfile(metrics_csv):
        return []
 
    with open(metrics_csv, "r", newline="") as f:
        reader = csv.DictReader(f)
        rows = list(reader)
    if not rows:
        return []
 
    x_name, x_values = infer_plot_x_axis(study_cfg, rows)
    if not x_name or x_values is None:
        print("[WARNING] Could not infer numeric x-axis for plots; skipping.")
        return []
 
    metric_keys = []
    param_keys = list((study_cfg.get("parameters") or {}).keys())
    for key in rows[0].keys():
        if key in {"case_id"}:
            continue
        if key in param_keys:
            continue
        metric_keys.append(key)
 
    out_format = plotting_cfg.get("output_format", "png")
    os.makedirs(plots_dir, exist_ok=True)
    generated = []
    for metric in metric_keys:
        y_values = []
        ok = True
        for row in rows:
            try:
                y_values.append(float(row[metric]))
            except Exception:
                ok = False
                break
        if not ok:
            continue
        plt.figure(figsize=(7.0, 4.2))
        plt.plot(x_values, y_values, marker="o", linewidth=1.5)
        plt.xlabel(x_name)
        plt.ylabel(metric)
        plt.title(f"{metric} vs {x_name}")
        plt.grid(True, alpha=0.3)
        out_path = os.path.join(plots_dir, f"{metric}_vs_{x_name.replace('/', '_')}.{out_format}")
        plt.tight_layout()
        plt.savefig(out_path, dpi=150)
        plt.close()
        generated.append(out_path)
    if generated:
        print(f"[SUCCESS] Generated {len(generated)} plot(s) in {os.path.relpath(plots_dir)}")
    return generated
 
 
def _command_to_string(command_tokens: list) -> str:
    """Render a command list as a shell-safe display string."""
    return " ".join(shlex.quote(str(tok)) for tok in command_tokens)
 
 
def _resolve_post_source_directory_preview(run_dir: str, monitor_cfg: dict, post_cfg: dict) -> str:
    """Resolve post source directory without side effects or stdout/stderr output."""
    solver_output_dir_rel = monitor_cfg.get('io', {}).get('directories', {}).get('output', 'results')
    solver_output_dir_abs = os.path.join(run_dir, solver_output_dir_rel)
    source_dir_template = post_cfg.get('source_data', {}).get('directory', '<solver_output_dir>')
    if source_dir_template == '<solver_output_dir>':
        return solver_output_dir_abs
    return os.path.abspath(os.path.join(run_dir, source_dir_template))
 
 
def build_run_dry_plan(args) -> dict:
    """Build a no-write execution plan for `run --dry-run`."""
    plan = {
        "mode": "dry-run",
        "created_at": datetime.now().isoformat(),
        "warnings": [],
        "inputs": {},
        "stages": {},
        "artifacts": [],
    }
 
    if args.dry_run and args.no_submit:
        plan["warnings"].append("--dry-run takes precedence over --no-submit; no files will be written.")
 
    cluster_mode = bool(getattr(args, "cluster", None))
    cluster_cfg = None
    cluster_path = None
    effective_num_procs = args.num_procs
    run_id = None
    run_dir = None
    solver_control_path = None
    loaded_case_cfg = None
    loaded_monitor_cfg = None
 
    if cluster_mode:
        cluster_path = os.path.abspath(args.cluster)
        cluster_cfg = read_yaml_file(cluster_path)
        validate_cluster_config(cluster_cfg, cluster_path)
        scheduler_type = str(cluster_cfg.get("scheduler", {}).get("type", "slurm")).lower()
        if args.scheduler and args.scheduler.lower() != scheduler_type:
            emit_structured_error(
                ERROR_CODE_CFG_INCONSISTENT_COMBO,
                key="scheduler.type",
                file_path=cluster_path,
                message=f"--scheduler={args.scheduler} does not match cluster.yml scheduler.type={scheduler_type}.",
            )
            sys.exit(1)
        if scheduler_type != "slurm":
            emit_structured_error(
                ERROR_CODE_CFG_INVALID_VALUE,
                key="scheduler.type",
                file_path=cluster_path,
                message=f"Unsupported scheduler '{scheduler_type}'. Only Slurm is supported in v1.",
            )
            sys.exit(1)
        cluster_tasks = get_cluster_total_tasks(cluster_cfg)
        if args.num_procs not in (1, cluster_tasks):
            emit_structured_error(
                ERROR_CODE_CFG_INCONSISTENT_COMBO,
                key="resources.ntasks_per_node",
                file_path=cluster_path,
                message=(
                    "--num-procs must be 1 (auto) or exactly nodes*ntasks_per_node "
                    f"({cluster_tasks}) in cluster mode."
                ),
            )
            sys.exit(1)
        effective_num_procs = cluster_tasks
        plan["launch_mode"] = "slurm"
        plan["inputs"]["cluster"] = cluster_path
    else:
        if getattr(args, "scheduler", None):
            fail_cli_usage("--scheduler requires --cluster in this version.")
        plan["launch_mode"] = "local"
 
    if args.solve:
        case_path = os.path.abspath(args.case)
        solver_path = os.path.abspath(args.solver)
        monitor_path = os.path.abspath(args.monitor)
        loaded_case_cfg = read_yaml_file(case_path)
        solver_cfg = read_yaml_file(solver_path)
        loaded_monitor_cfg = read_yaml_file(monitor_path)
        validate_solver_configs(loaded_case_cfg, solver_cfg, loaded_monitor_cfg, case_path, solver_path, monitor_path)
 
        case_name = os.path.splitext(os.path.basename(case_path))[0]
        timestamp = datetime.now().strftime("%Y%m%d-%H%M%S")
        run_id = f"{case_name}_{timestamp}"
        run_dir = os.path.abspath(os.path.join("runs", run_id))
 
        config_dir = os.path.join(run_dir, "config")
        scheduler_dir = os.path.join(run_dir, "scheduler")
        logs_dir = os.path.join(run_dir, "logs")
        solver_control_path = os.path.join(config_dir, f"{run_id}.control")
        whitelist_path = os.path.join(config_dir, "whitelist.run")
        profile_path = os.path.join(config_dir, "profile.run")
 
        plan["run_id_preview"] = run_id
        plan["run_dir_preview"] = run_dir
        plan["inputs"].update({"case": case_path, "solver": solver_path, "monitor": monitor_path})
        plan["artifacts"].extend(
            [
                run_dir,
                config_dir,
                logs_dir,
                os.path.join(run_dir, "results"),
                scheduler_dir,
                os.path.join(config_dir, "case.yml"),
                os.path.join(config_dir, "solver.yml"),
                os.path.join(config_dir, "monitor.yml"),
                whitelist_path,
                profile_path,
                solver_control_path,
                os.path.join(run_dir, "manifest.json"),
            ]
        )
        if cluster_mode:
            plan["artifacts"].append(os.path.join(config_dir, "cluster.yml"))
            plan["artifacts"].append(os.path.join(scheduler_dir, "submission.json"))
 
        solver_exe = os.path.join(BIN_DIR, "picsolver")
        solver_args = ["-control_file", solver_control_path]
        if cluster_mode:
            solver_script = os.path.join(scheduler_dir, "solver.sbatch")
            solver_cmd = build_cluster_launch_command(cluster_cfg, solver_exe, solver_args)
            plan["artifacts"].append(solver_script)
            plan["stages"]["solve"] = {
                "mode": "slurm",
                "script": solver_script,
                "launch_command": solver_cmd,
                "launch_command_string": _command_to_string(solver_cmd),
            }
        else:
            solver_cmd = [solver_exe] + solver_args
            if effective_num_procs > 1:
                solver_cmd = ["mpiexec", "-n", str(effective_num_procs)] + solver_cmd
            plan["stages"]["solve"] = {
                "mode": "local",
                "launch_command": solver_cmd,
                "launch_command_string": _command_to_string(solver_cmd),
            }
 
    if args.post_process:
        post_path = os.path.abspath(args.post)
        plan["inputs"]["post"] = post_path
        post_cfg = read_yaml_file(post_path)
        validate_post_config(post_cfg, post_path)
 
        if args.run_dir:
            run_dir = os.path.abspath(args.run_dir)
            if not os.path.isdir(run_dir):
                emit_structured_error(
                    ERROR_CODE_CFG_FILE_NOT_FOUND,
                    key="run-dir",
                    file_path=run_dir,
                    message="Specified run directory not found.",
                )
                sys.exit(1)
            run_id = os.path.basename(run_dir)
        elif not args.solve:
            fail_cli_usage("--post-process requires --run-dir when not used with --solve.")
 
        if args.run_dir:
            config_dir = os.path.join(run_dir, "config")
            case_path, monitor_path, solver_control_path = auto_identify_run_inputs(config_dir)
            if not all([case_path, monitor_path, solver_control_path]):
                emit_structured_error(
                    ERROR_CODE_CFG_MISSING_KEY,
                    key="run_dir.config",
                    file_path=config_dir,
                    message=(
                        "Could not auto-identify required run inputs "
                        "(case.yml/monitor.yml/*.control) in run config directory."
                    ),
                )
                sys.exit(1)
            loaded_case_cfg = read_yaml_file(case_path)
            loaded_monitor_cfg = read_yaml_file(monitor_path)
        else:
            config_dir = os.path.join(run_dir, "config")
            case_path = os.path.join(config_dir, "case.yml")
            monitor_path = os.path.join(config_dir, "monitor.yml")
            if solver_control_path is None:
                solver_control_path = os.path.join(config_dir, f"{run_id}.control")
 
        resolved_source_dir = _resolve_post_source_directory_preview(run_dir, loaded_monitor_cfg, post_cfg)
        post_recipe_path = os.path.join(config_dir, "post.run")
        output_dir_rel = post_cfg.get("io", {}).get("output_directory")
        output_prefix = post_cfg.get("io", {}).get("output_filename_prefix")
        if not output_dir_rel or not output_prefix:
            emit_structured_error(
                ERROR_CODE_CFG_MISSING_KEY,
                key="io.output_directory/io.output_filename_prefix",
                file_path=post_path,
                message="Missing required post IO keys.",
            )
            sys.exit(1)
        output_dir_abs = os.path.abspath(os.path.join(run_dir, output_dir_rel))
        post_exe = os.path.join(BIN_DIR, "postprocessor")
        post_args = ["-control_file", solver_control_path, "-postprocessing_config_file", post_recipe_path]
        plan["artifacts"].extend([post_recipe_path, output_dir_abs])
 
        if cluster_mode:
            scheduler_dir = os.path.join(run_dir, "scheduler")
            post_script = os.path.join(scheduler_dir, "post.sbatch")
            post_cmd = build_cluster_launch_command(cluster_cfg, post_exe, post_args)
            plan["artifacts"].append(post_script)
            plan["stages"]["post-process"] = {
                "mode": "slurm",
                "script": post_script,
                "source_data_directory": resolved_source_dir,
                "launch_command": post_cmd,
                "launch_command_string": _command_to_string(post_cmd),
            }
        else:
            post_cmd = [post_exe] + post_args
            if effective_num_procs > 1:
                post_cmd = ["mpiexec", "-n", str(effective_num_procs)] + post_cmd
            plan["stages"]["post-process"] = {
                "mode": "local",
                "source_data_directory": resolved_source_dir,
                "launch_command": post_cmd,
                "launch_command_string": _command_to_string(post_cmd),
            }
 
    # Preserve insertion order while removing duplicates.
    deduped = []
    seen = set()
    for item in plan["artifacts"]:
        if item not in seen:
            seen.add(item)
            deduped.append(item)
    plan["artifacts"] = deduped
    if run_id and "run_id_preview" not in plan:
        plan["run_id_preview"] = run_id
    if run_dir and "run_dir_preview" not in plan:
        plan["run_dir_preview"] = run_dir
    plan["num_procs_effective"] = effective_num_procs
    return plan
 
 
def render_run_dry_plan(plan: dict, output_format: str = "text"):
    """Render dry-run plan in human or JSON format."""
    if output_format == "json":
        print(json.dumps(plan, indent=2, sort_keys=True))
        return
 
    print("\n" + "=" * 60)
    print("                      DRY-RUN PLAN")
    print("=" * 60)
    print(f"  Launch mode    : {plan.get('launch_mode')}")
    print(f"  Created at     : {plan.get('created_at')}")
    if plan.get("run_id_preview"):
        print(f"  Run ID preview : {plan.get('run_id_preview')}")
    if plan.get("run_dir_preview"):
        print(f"  Run dir preview: {plan.get('run_dir_preview')}")
    print(f"  MPI processes  : {plan.get('num_procs_effective')}")
    if plan.get("warnings"):
        print("  Warnings       :")
        for warning in plan["warnings"]:
            print(f"    - {warning}")
 
    if plan.get("inputs"):
        print("\n  Inputs:")
        for key, value in plan["inputs"].items():
            print(f"    - {key}: {value}")
 
    if plan.get("stages"):
        print("\n  Planned stage commands:")
        for stage, details in plan["stages"].items():
            print(f"    - {stage} ({details.get('mode')}):")
            print(f"      {details.get('launch_command_string')}")
 
    print("\n  Planned artifacts (no files created in dry-run):")
    for artifact in plan.get("artifacts", []):
        print(f"    - {artifact}")
    print("=" * 60)
 
 
def validate_workflow(args):
    """Implements `pic.flow validate` without launching solver/post workflows."""
    checked = []
    solver_group_selected = any([args.case, args.solver, args.monitor])
    any_group_selected = solver_group_selected or any([args.post, args.cluster, args.study])
 
    if not any_group_selected:
        fail_cli_usage(
            "validate requires at least one config group. Provide solver trio and/or --post/--cluster/--study.",
            hint="Example: pic.flow validate --case case.yml --solver solver.yml --monitor monitor.yml --post post.yml",
        )
 
    if solver_group_selected and not all([args.case, args.solver, args.monitor]):
        fail_cli_usage("When solver validation is requested, --case, --solver, and --monitor are all required.")
 
    if solver_group_selected:
        case_path = os.path.abspath(args.case)
        solver_path = os.path.abspath(args.solver)
        monitor_path = os.path.abspath(args.monitor)
        case_cfg = read_yaml_file(case_path)
        solver_cfg = read_yaml_file(solver_path)
        monitor_cfg = read_yaml_file(monitor_path)
        validate_solver_configs(case_cfg, solver_cfg, monitor_cfg, case_path, solver_path, monitor_path)
        checked.extend([case_path, solver_path, monitor_path])
 
    post_cfg = None
    if args.post:
        post_path = os.path.abspath(args.post)
        post_cfg = read_yaml_file(post_path)
        validate_post_config(post_cfg, post_path)
        checked.append(post_path)
 
    cluster_cfg = None
    if args.cluster:
        cluster_path = os.path.abspath(args.cluster)
        cluster_cfg = read_yaml_file(cluster_path)
        validate_cluster_config(cluster_cfg, cluster_path)
        checked.append(cluster_path)
 
    study_cfg = None
    if args.study:
        study_path = os.path.abspath(args.study)
        study_cfg = read_yaml_file(study_path)
        validate_study_config(study_cfg, study_path)
        checked.append(study_path)
 
    if args.strict and post_cfg is not None:
        post_path = os.path.abspath(args.post)
        source_dir = post_cfg.get("source_data", {}).get("directory")
        if source_dir and source_dir != "<solver_output_dir>":
            resolved = resolve_path(post_path, source_dir)
            if not os.path.isdir(resolved):
                emit_structured_error(
                    ERROR_CODE_CFG_FILE_NOT_FOUND,
                    key="source_data.directory",
                    file_path=post_path,
                    message=f"strict mode: source_data.directory resolves to missing directory '{resolved}'.",
                )
                sys.exit(1)
 
    if args.strict and study_cfg is not None:
        study_path = os.path.abspath(args.study)
        base_cfgs = study_cfg.get("base_configs", {})
        if isinstance(base_cfgs, dict):
            base_case_path = resolve_path(study_path, base_cfgs.get("case"))
            base_solver_path = resolve_path(study_path, base_cfgs.get("solver"))
            base_monitor_path = resolve_path(study_path, base_cfgs.get("monitor"))
            base_post_path = resolve_path(study_path, base_cfgs.get("post"))
            if all([base_case_path, base_solver_path, base_monitor_path]):
                validate_solver_configs(
                    read_yaml_file(base_case_path),
                    read_yaml_file(base_solver_path),
                    read_yaml_file(base_monitor_path),
                    base_case_path,
                    base_solver_path,
                    base_monitor_path,
                )
            if base_post_path:
                validate_post_config(read_yaml_file(base_post_path), base_post_path)
 
    print(f"[SUCCESS] Validation completed for {len(checked)} file(s).")
    for path in checked:
        print(f"  - {path}")
 
def run_workflow(args):
    """Main orchestrator for the 'run' command (local and Slurm modes)."""
    if getattr(args, "dry_run", False):
        plan = build_run_dry_plan(args)
        render_run_dry_plan(plan, output_format=getattr(args, "output_format", "text"))
        return
 
    run_dir = None
    run_id = None
    output_dir_abs = None
    workflow_start = time.time()
    stages_completed = []
    configs = None
    submission_meta = {"launch_mode": "local", "stages": {}}
 
    cluster_mode = bool(getattr(args, "cluster", None))
    cluster_cfg = None
    cluster_path = None
    effective_num_procs = args.num_procs
 
    if cluster_mode:
        cluster_path = os.path.abspath(args.cluster)
        cluster_cfg = read_yaml_file(cluster_path)
        validate_cluster_config(cluster_cfg, cluster_path)
        scheduler_type = str(cluster_cfg.get("scheduler", {}).get("type", "slurm")).lower()
        if args.scheduler and args.scheduler.lower() != scheduler_type:
            print(
                f"[FATAL] --scheduler={args.scheduler} does not match cluster.yml scheduler.type={scheduler_type}.",
                file=sys.stderr
            )
            sys.exit(1)
        if scheduler_type != "slurm":
            print(f"[FATAL] Unsupported scheduler '{scheduler_type}'. Only Slurm is supported in v1.", file=sys.stderr)
            sys.exit(1)
        cluster_tasks = get_cluster_total_tasks(cluster_cfg)
        if args.num_procs not in (1, cluster_tasks):
            print(
                f"[FATAL] In cluster mode, --num-procs must be 1 (auto) or exactly nodes*ntasks_per_node ({cluster_tasks}).",
                file=sys.stderr
            )
            sys.exit(1)
        effective_num_procs = cluster_tasks
        submission_meta["launch_mode"] = "slurm"
        submission_meta["cluster_config"] = cluster_path
        submission_meta["no_submit"] = bool(args.no_submit)
        print(f"[INFO] Cluster mode enabled (Slurm). Using {effective_num_procs} MPI tasks from cluster.yml.")
    elif getattr(args, "scheduler", None):
        print("[FATAL] --scheduler requires --cluster in this version.", file=sys.stderr)
        sys.exit(1)
 
    # --- Stage 1: Solver (if requested) ---
    if args.solve:
        configs = {
            'case': read_yaml_file(args.case), 'case_path': os.path.abspath(args.case),
            'solver': read_yaml_file(args.solver), 'solver_path': os.path.abspath(args.solver),
            'monitor': read_yaml_file(args.monitor), 'monitor_path': os.path.abspath(args.monitor)
        }
 
        print("\n[INFO] Validating configuration files...")
        validate_solver_configs(
            configs['case'], configs['solver'], configs['monitor'],
            args.case, args.solver, args.monitor
        )
        print("[SUCCESS] All configuration files passed validation.\n")
 
        case_name = os.path.splitext(os.path.basename(args.case))[0]
        timestamp = datetime.now().strftime("%Y%m%d-%H%M%S")
        run_id = f"{case_name}_{timestamp}"
        run_dir = os.path.abspath(os.path.join("runs", run_id))
 
        config_dir = os.path.join(run_dir, "config")
        for d in [config_dir, "logs", "results", os.path.join(run_dir, "scheduler")]:
            os.makedirs(d, exist_ok=True)
        print(f"[INFO] Created new self-contained run directory: {os.path.relpath(run_dir)}")
 
        shutil.copy(args.case, os.path.join(config_dir, "case.yml"))
        shutil.copy(args.solver, os.path.join(config_dir, "solver.yml"))
        shutil.copy(args.monitor, os.path.join(config_dir, "monitor.yml"))
        if cluster_mode:
            shutil.copy(cluster_path, os.path.join(config_dir, "cluster.yml"))
 
        print("\n" + "="*25 + " SOLVER STAGE " + "="*25)
        source_files = {'Case': args.case, 'Solver': args.solver, 'Monitor': args.monitor}
        print("[INFO] Generating monitoring files (whitelist.run, profile.run)...")
        whitelist_path = generate_simple_list_file(
            run_dir, run_id, configs['monitor'], 'logging', 'enabled_functions', 'whitelist.run', source_files
        )
        profile_path = generate_simple_list_file(
            run_dir, run_id, configs['monitor'], 'profiling', 'critical_functions', 'profile.run', source_files
        )
 
        monitor_files = {'whitelist': whitelist_path, 'profile': profile_path}
        control_file = generate_solver_control_file(run_dir, run_id, configs, effective_num_procs, monitor_files)
 
        solver_exe = os.path.join(BIN_DIR, "picsolver")
        solver_args = ["-control_file", control_file]
        if cluster_mode:
            scheduler_dir = os.path.join(run_dir, "scheduler")
            solver_script = os.path.join(scheduler_dir, "solver.sbatch")
            solver_log = os.path.join(run_dir, "logs", "solver_%j.out")
            solver_err = os.path.join(run_dir, "logs", "solver_%j.err")
            solver_cmd = build_cluster_launch_command(cluster_cfg, solver_exe, solver_args)
            render_slurm_script(
                solver_script,
                f"{run_id}_solve",
                cluster_cfg,
                solver_cmd,
                run_dir,
                solver_log,
                solver_err,
                env_vars={"LOG_LEVEL": configs['monitor'].get('logging', {}).get('verbosity', 'INFO').upper()},
            )
            submission_meta["stages"]["solve"] = {"script": solver_script, "submitted": False}
            print(f"[SUCCESS] Generated solver Slurm script: {os.path.relpath(solver_script)}")
            if not args.no_submit:
                submit_info = submit_sbatch(solver_script)
                submission_meta["stages"]["solve"].update(submit_info)
                submission_meta["stages"]["solve"]["submitted"] = True
                print(f"[SUCCESS] Submitted solver job: {submit_info['job_id']}")
            stages_completed.append('solve')
        else:
            command = [solver_exe] + solver_args
            if effective_num_procs > 1:
                command = ["mpiexec", "-n", str(effective_num_procs)] + command
            execute_command(command, run_dir, f"{run_id}_solver.log", configs['monitor'])
            stages_completed.append('solve')
 
    # --- Stage 2: Post-Processing (if requested) ---
    if args.post_process:
        if args.run_dir:
            run_dir = os.path.abspath(args.run_dir)
            if not os.path.isdir(run_dir):
                print(f"[FATAL] Specified run directory not found: {run_dir}", file=sys.stderr)
                sys.exit(1)
            print(f"[INFO] Operating on existing run directory: {os.path.relpath(run_dir)}")
            run_id = os.path.basename(run_dir)
        elif not args.solve:
            print("[FATAL] --post-process requires --run-dir when not used with --solve.", file=sys.stderr)
            sys.exit(1)
 
        print("\n" + "="*20 + " POST-PROCESSING STAGE " + "="*20)
        config_dir = os.path.join(run_dir, "config")
        case_path, monitor_path, solver_control_path = auto_identify_run_inputs(config_dir)
 
        if not all([case_path, monitor_path, solver_control_path]):
            print(f"[FATAL] Could not automatically identify required config files in {config_dir}", file=sys.stderr)
            if not case_path:
                print("         - No 'case' file found (expected 'models' + 'boundary_conditions').", file=sys.stderr)
            if not monitor_path:
                print("         - No 'monitor' file found (expected 'io' + 'logging').", file=sys.stderr)
            if not solver_control_path:
                print("         - No '.control' file found.", file=sys.stderr)
            sys.exit(1)
 
        print(f"[INFO] Auto-identified Case file:    {os.path.basename(case_path)}")
        print(f"[INFO] Auto-identified Monitor file: {os.path.basename(monitor_path)}")
 
        case_cfg = read_yaml_file(case_path)
        monitor_cfg = read_yaml_file(monitor_path)
        post_cfg = read_yaml_file(args.post)
 
        print("[INFO] Validating post-processing configuration...")
        validate_post_config(post_cfg, args.post)
        print("[SUCCESS] Post-processing configuration passed validation.\n")
 
        strict_source_check = not (cluster_mode and args.solve)
        resolved_source_dir = resolve_post_source_directory(run_dir, monitor_cfg, post_cfg, strict=strict_source_check)
        if 'source_data' not in post_cfg:
            post_cfg['source_data'] = {}
        post_cfg['source_data']['directory'] = resolved_source_dir
 
        post_io_cfg = post_cfg.get('io', {})
        try:
            output_dir_rel = post_io_cfg['output_directory']
            output_prefix = post_io_cfg['output_filename_prefix']
        except KeyError as e:
            print(f"[FATAL] Missing required key '{e.args[0]}' in the 'io' section of {args.post}", file=sys.stderr)
            sys.exit(1)
 
        output_dir_abs = os.path.abspath(os.path.join(run_dir, output_dir_rel))
        os.makedirs(output_dir_abs, exist_ok=True)
        print(f"[INFO] Post-processor output directory: {os.path.relpath(output_dir_abs)}")
 
        source_files_post = {'Case': case_path, 'Post-Profile': args.post}
        post_recipe_file = generate_post_recipe_file(run_dir, run_id, post_cfg, source_files_post)
 
        post_exe = os.path.join(BIN_DIR, "postprocessor")
        post_args = ["-control_file", solver_control_path, "-postprocessing_config_file", post_recipe_file]
        if cluster_mode:
            scheduler_dir = os.path.join(run_dir, "scheduler")
            os.makedirs(scheduler_dir, exist_ok=True)
            post_script = os.path.join(scheduler_dir, "post.sbatch")
            post_log = os.path.join(run_dir, "logs", "post_%j.out")
            post_err = os.path.join(run_dir, "logs", "post_%j.err")
            post_cmd = build_cluster_launch_command(cluster_cfg, post_exe, post_args)
            render_slurm_script(
                post_script,
                f"{run_id}_post",
                cluster_cfg,
                post_cmd,
                run_dir,
                post_log,
                post_err,
                env_vars={"LOG_LEVEL": monitor_cfg.get('logging', {}).get('verbosity', 'INFO').upper()},
            )
            submission_meta["stages"]["post-process"] = {"script": post_script, "submitted": False}
            print(f"[SUCCESS] Generated post Slurm script: {os.path.relpath(post_script)}")
 
            if not args.no_submit:
                dependency_job = None
                if args.solve:
                    dependency_job = submission_meta.get("stages", {}).get("solve", {}).get("job_id")
                submit_info = submit_sbatch(post_script, dependency=dependency_job)
                submission_meta["stages"]["post-process"].update(submit_info)
                submission_meta["stages"]["post-process"]["submitted"] = True
                if dependency_job:
                    submission_meta["stages"]["post-process"]["dependency"] = f"afterok:{dependency_job}"
                print(f"[SUCCESS] Submitted post job: {submit_info['job_id']}")
            stages_completed.append('post-process')
        else:
            command = [post_exe] + post_args
            if effective_num_procs > 1:
                command = ["mpiexec", "-n", str(effective_num_procs)] + command
            execute_command(command, run_dir, f"{run_id}_{output_prefix}.log", monitor_cfg)
            stages_completed.append('post-process')
 
    if run_dir:
        manifest = {
            "run_id": run_id,
            "created_at": datetime.now().isoformat(),
            "launch_mode": "slurm" if cluster_mode else "local",
            "git_commit": get_git_commit(),
            "num_procs": effective_num_procs,
            "stages_requested": {"solve": bool(args.solve), "post_process": bool(args.post_process)},
            "stages_completed_or_submitted": stages_completed,
            "inputs": {},
        }
        if args.solve:
            manifest["inputs"]["case"] = os.path.abspath(args.case)
            manifest["inputs"]["solver"] = os.path.abspath(args.solver)
            manifest["inputs"]["monitor"] = os.path.abspath(args.monitor)
        if args.post_process:
            manifest["inputs"]["post"] = os.path.abspath(args.post)
        if cluster_mode:
            manifest["inputs"]["cluster"] = cluster_path
            write_json_file(os.path.join(run_dir, "scheduler", "submission.json"), submission_meta)
        write_json_file(os.path.join(run_dir, "manifest.json"), manifest)
 
    if stages_completed:
        elapsed = time.time() - workflow_start
        mins, secs = divmod(int(elapsed), 60)
        hrs, mins = divmod(mins, 60)
        if hrs > 0:
            time_str = f"{hrs}h {mins}m {secs}s"
        elif mins > 0:
            time_str = f"{mins}m {secs}s"
        else:
            time_str = f"{secs}s"
 
        print("\n" + "=" * 60)
        print("                       RUN SUMMARY")
        print("=" * 60)
        print(f"  Run ID         : {run_id}")
        print(f"  Run directory  : {os.path.relpath(run_dir)}")
        print(f"  Wall-clock     : {time_str}")
        print(f"  Stages         : {', '.join(stages_completed)}")
        print(f"  Launch mode    : {'slurm' if cluster_mode else 'local'}")
        print(f"  MPI processes  : {effective_num_procs}")
        if args.solve and configs:
            total_steps = configs['case'].get('run_control', {}).get('total_steps', '?')
            result_dir = os.path.join(run_dir, configs['monitor'].get('io', {}).get('directories', {}).get('output', 'results'))
            print(f"  Steps run      : {total_steps}")
            print(f"  Solver output  : {os.path.relpath(result_dir)}")
        if 'post-process' in stages_completed and output_dir_abs:
            print(f"  Post output    : {os.path.relpath(output_dir_abs)}")
        print(f"  Logs           : {os.path.relpath(os.path.join(run_dir, 'logs'))}")
        if cluster_mode:
            submission_file = os.path.join(run_dir, "scheduler", "submission.json")
            print(f"  Submission meta: {os.path.relpath(submission_file)}")
        print("=" * 60)
 
def sweep_workflow(args):
    """Study/sweep orchestration using Slurm job arrays."""
    study_path = os.path.abspath(args.study)
    cluster_path = os.path.abspath(args.cluster)
 
    study_cfg = read_yaml_file(study_path)
    cluster_cfg = read_yaml_file(cluster_path)
    validate_study_config(study_cfg, study_path)
    validate_cluster_config(cluster_cfg, cluster_path)
 
    study_name = os.path.splitext(os.path.basename(study_path))[0]
    timestamp = datetime.now().strftime("%Y%m%d-%H%M%S")
    study_id = f"{study_name}_{timestamp}"
    study_dir = os.path.abspath(os.path.join("studies", study_id))
    cases_dir = os.path.join(study_dir, "cases")
    scheduler_dir = os.path.join(study_dir, "scheduler")
    results_dir = os.path.join(study_dir, "results")
    logs_dir = os.path.join(study_dir, "logs")
    for path in [cases_dir, scheduler_dir, results_dir, logs_dir]:
        os.makedirs(path, exist_ok=True)
 
    print(f"[INFO] Creating study directory: {os.path.relpath(study_dir)}")
    shutil.copy(study_path, os.path.join(study_dir, "study.yml"))
    shutil.copy(cluster_path, os.path.join(study_dir, "cluster.yml"))
 
    base_cfgs = study_cfg["base_configs"]
    base_paths = {k: resolve_path(study_path, v) for k, v in base_cfgs.items()}
    base_case = read_yaml_file(base_paths["case"])
    base_solver = read_yaml_file(base_paths["solver"])
    base_monitor = read_yaml_file(base_paths["monitor"])
    base_post = read_yaml_file(base_paths["post"])
    validate_solver_configs(base_case, base_solver, base_monitor, base_paths["case"], base_paths["solver"], base_paths["monitor"])
    validate_post_config(base_post, base_paths["post"])
 
    combinations = expand_parameter_matrix(study_cfg["parameters"])
    if not combinations:
        print("[FATAL] Study parameter matrix expanded to zero cases.", file=sys.stderr)
        sys.exit(1)
    print(f"[INFO] Expanded sweep matrix to {len(combinations)} case(s).")
 
    cluster_tasks = get_cluster_total_tasks(cluster_cfg)
    case_entries = []
    case_index_file = os.path.join(scheduler_dir, "case_index.tsv")
 
    for idx, combo in enumerate(combinations):
        case_id = f"case_{idx:04d}"
        run_dir = os.path.join(cases_dir, case_id)
        config_dir = os.path.join(run_dir, "config")
        os.makedirs(config_dir, exist_ok=True)
        os.makedirs(os.path.join(run_dir, "logs"), exist_ok=True)
        os.makedirs(os.path.join(run_dir, "results"), exist_ok=True)
 
        case_cfg = copy.deepcopy(base_case)
        solver_cfg = copy.deepcopy(base_solver)
        monitor_cfg = copy.deepcopy(base_monitor)
        post_cfg = copy.deepcopy(base_post)
        target_map = {"case": case_cfg, "solver": solver_cfg, "monitor": monitor_cfg, "post": post_cfg}
        for full_key, value in combo.items():
            root, nested = full_key.split(".", 1)
            _deep_set(target_map[root], nested, value)
 
        # Preserve file-based/grid-gen workflows when study cases are materialized
        # into new directories by rewriting external paths as absolute.
        absolutize_case_external_paths(case_cfg, base_paths["case"])
 
        case_path = os.path.join(config_dir, "case.yml")
        solver_path = os.path.join(config_dir, "solver.yml")
        monitor_path = os.path.join(config_dir, "monitor.yml")
        post_path = os.path.join(config_dir, "post.yml")
        write_yaml_file(case_path, case_cfg)
        write_yaml_file(solver_path, solver_cfg)
        write_yaml_file(monitor_path, monitor_cfg)
        write_yaml_file(post_path, post_cfg)
 
        validate_solver_configs(case_cfg, solver_cfg, monitor_cfg, case_path, solver_path, monitor_path)
        validate_post_config(post_cfg, post_path)
 
        source_files = {'Case': case_path, 'Solver': solver_path, 'Monitor': monitor_path}
        whitelist_path = generate_simple_list_file(run_dir, case_id, monitor_cfg, 'logging', 'enabled_functions', 'whitelist.run', source_files)
        profile_path = generate_simple_list_file(run_dir, case_id, monitor_cfg, 'profiling', 'critical_functions', 'profile.run', source_files)
        monitor_files = {"whitelist": whitelist_path, "profile": profile_path}
        configs = {
            "case": case_cfg, "case_path": case_path,
            "solver": solver_cfg, "solver_path": solver_path,
            "monitor": monitor_cfg, "monitor_path": monitor_path
        }
        control_file = generate_solver_control_file(run_dir, case_id, configs, cluster_tasks, monitor_files)
 
        source_dir = resolve_post_source_directory(run_dir, monitor_cfg, post_cfg, strict=False)
        if 'source_data' not in post_cfg:
            post_cfg['source_data'] = {}
        post_cfg['source_data']['directory'] = source_dir
        output_prefix = post_cfg.get("io", {}).get("output_filename_prefix", "post")
        post_recipe = generate_post_recipe_file(run_dir, case_id, post_cfg, {'Case': case_path, 'Post-Profile': post_path})
 
        case_entries.append({
            "index": idx,
            "case_id": case_id,
            "run_dir": os.path.abspath(run_dir),
            "control_file": control_file,
            "post_recipe_file": post_recipe,
            "log_level": str(monitor_cfg.get("logging", {}).get("verbosity", "INFO")).upper(),
            "post_prefix": output_prefix,
            "parameters": combo,
        })
 
    with open(case_index_file, "w") as f:
        for entry in case_entries:
            f.write(
                "\t".join(
                    [
                        str(entry["index"]),
                        entry["case_id"],
                        entry["run_dir"],
                        entry["control_file"],
                        entry["post_recipe_file"],
                        entry["log_level"],
                        entry["post_prefix"],
                    ]
                ) + "\n"
            )
    print(f"[SUCCESS] Wrote sweep case index: {os.path.relpath(case_index_file)}")
 
    max_idx = len(case_entries) - 1
    max_conc = study_cfg.get("execution", {}).get("max_concurrent_array_tasks")
    array_spec = f"0-{max_idx}"
    if max_conc:
        array_spec = f"{array_spec}%{max_conc}"
 
    solver_exe = os.path.join(BIN_DIR, "picsolver")
    post_exe = os.path.join(BIN_DIR, "postprocessor")
    solver_array_script = os.path.join(scheduler_dir, "solver_array.sbatch")
    post_array_script = os.path.join(scheduler_dir, "post_array.sbatch")
    render_slurm_array_stage_script(
        solver_array_script,
        f"{study_id}_solve",
        cluster_cfg,
        array_spec,
        case_index_file,
        "solve",
        solver_exe,
        post_exe,
        os.path.join(logs_dir, "solver_%A_%a.out"),
        os.path.join(logs_dir, "solver_%A_%a.err")
    )
    render_slurm_array_stage_script(
        post_array_script,
        f"{study_id}_post",
        cluster_cfg,
        array_spec,
        case_index_file,
        "post",
        solver_exe,
        post_exe,
        os.path.join(logs_dir, "post_%A_%a.out"),
        os.path.join(logs_dir, "post_%A_%a.err")
    )
    print(f"[SUCCESS] Generated Slurm array scripts in {os.path.relpath(scheduler_dir)}")
 
    submission = {
        "launch_mode": "slurm",
        "study_id": study_id,
        "solver_array": {"script": solver_array_script, "submitted": False},
        "post_array": {"script": post_array_script, "submitted": False},
        "no_submit": bool(args.no_submit),
    }
    if not args.no_submit:
        solver_submit = submit_sbatch(solver_array_script)
        submission["solver_array"].update(solver_submit)
        submission["solver_array"]["submitted"] = True
        post_submit = submit_sbatch(post_array_script, dependency=solver_submit["job_id"])
        submission["post_array"].update(post_submit)
        submission["post_array"]["submitted"] = True
        submission["post_array"]["dependency"] = f"afterok:{solver_submit['job_id']}"
        print(f"[SUCCESS] Submitted solver array job: {solver_submit['job_id']}")
        print(f"[SUCCESS] Submitted post array job:   {post_submit['job_id']}")
 
    metrics_csv = aggregate_study_metrics(study_cfg, case_entries, results_dir)
    plots = generate_study_plots(study_cfg, metrics_csv, os.path.join(results_dir, "plots"))
 
    summary = {
        "study_id": study_id,
        "created_at": datetime.now().isoformat(),
        "git_commit": get_git_commit(),
        "study_type": study_cfg.get("study_type"),
        "num_cases": len(case_entries),
        "paths": {
            "study_dir": study_dir,
            "case_index": case_index_file,
            "solver_array_script": solver_array_script,
            "post_array_script": post_array_script,
            "metrics_table": metrics_csv,
            "plots_dir": os.path.join(results_dir, "plots"),
        },
        "submission": submission,
    }
    write_json_file(os.path.join(scheduler_dir, "submission.json"), submission)
    write_json_file(os.path.join(study_dir, "study_manifest.json"), summary)
    write_json_file(os.path.join(results_dir, "summary.json"), {"study_id": study_id, "metrics_csv": metrics_csv, "plots": plots})
 
    print("\n" + "=" * 60)
    print("                     STUDY SUMMARY")
    print("=" * 60)
    print(f"  Study ID        : {study_id}")
    print(f"  Study directory : {os.path.relpath(study_dir)}")
    print(f"  Cases generated : {len(case_entries)}")
    print(f"  Array spec      : {array_spec}")
    print(f"  Solver script   : {os.path.relpath(solver_array_script)}")
    print(f"  Post script     : {os.path.relpath(post_array_script)}")
    if metrics_csv:
        print(f"  Metrics table   : {os.path.relpath(metrics_csv)}")
    if plots:
        print(f"  Plots           : {os.path.relpath(os.path.join(results_dir, 'plots'))}")
    print("=" * 60)
 
 
def init_case(args):
    """!
    @brief Implements the 'init' command.
    @details Creates a new case study directory by copying a template. It can then
             either create relative symbolic links to the project's executables
             (default) or create a full copy of them for a self-contained study.
    @param[in] args The command-line arguments parsed by argparse.
    """
    template_path = os.path.join(PROJECT_ROOT, "examples", args.template_name)
    # The destination path is relative to the current working directory.
    dest_path = os.path.abspath(os.path.join(os.getcwd(), args.dest_name if args.dest_name else args.template_name))
 
    if not os.path.isdir(template_path):
        print(f"[FATAL] Case template '{args.template_name}' not found at '{template_path}'", file=sys.stderr)
        sys.exit(1)
    if os.path.exists(dest_path):
        print(f"[FATAL] Destination directory '{dest_path}' already exists.", file=sys.stderr)
        sys.exit(1)
 
    print(f"[INFO] Initializing new case '{os.path.basename(dest_path)}' from template '{args.template_name}'...")
    
    shutil.copytree(template_path, dest_path)
    print(f"[SUCCESS] Copied template files to: {dest_path}")
    
    if args.copy_binaries:
        print("[INFO] Copying project binaries for a self-contained case...")
    else:
        print("[INFO] Creating symbolic links to project binaries...")
        
    main_bin_dir = os.path.join(PROJECT_ROOT, "bin")
    
    try:
        # We only want to copy/link the actual executables, not the orchestrator script itself.
        binaries = [f for f in os.listdir(main_bin_dir) if f != 'pic.flow' and os.path.isfile(os.path.join(main_bin_dir, f))]
        if not binaries:
            print("[WARNING] Main project bin/ directory contains no executables. Nothing to link or copy.", file=sys.stderr)
            return
 
        for binary_name in binaries:
            source_path_abs = os.path.abspath(os.path.join(main_bin_dir, binary_name))
            dest_file_path = os.path.join(dest_path, binary_name)
            
            if args.copy_binaries:
                # Use copy2 to preserve permissions and metadata
                shutil.copy2(source_path_abs, dest_file_path)
                print(f"  - Copied '{binary_name}'")
            else:
                relative_source_path = os.path.relpath(source_path_abs, start=dest_path)
                os.symlink(relative_source_path, dest_file_path)
                print(f"  - Linked '{binary_name}'")
            
        print("[SUCCESS] Binaries are now available in your study directory.")
        print("          You can now 'cd' into your study and run commands locally (e.g., './picsolver ...').")
 
    except Exception as e:
        print(f"[ERROR] Failed to process binaries: {e}", file=sys.stderr)
        print("        Your case files were copied, but you will need to run commands from the project root.", file=sys.stderr)
 
def build_project(args):
    """!
    @brief Implements the 'build' command.
    @details Executes the top-level build.sh script, passing through any
             additional arguments directly to 'make'. This allows for building,
             cleaning, and other Makefile targets via the orchestrator.
    @param[in] args The command-line arguments parsed by argparse.
    """
 
    print("\n" + "="*27 + " BUILD STAGE " + "="*27)
    build_script_path = os.path.join(PROJECT_ROOT, "build.sh")
 
    if not os.path.isfile(build_script_path):
        print(f"[FATAL] Build script not found at expected location: {build_script_path}", file=sys.stderr)
        print("        Please ensure 'build.sh' exists in the project root directory.", file=sys.stderr)
        sys.exit(1)
 
    # The command is the script itself, plus any passthrough arguments for make.
   # command = [build_script_path] + args.make_args
    command = ['/bin/bash', build_script_path] + args.make_args
    # For the build process, we don't have a monitor.yml, so we pass an empty
    # dict to execute_command. The command should be run in the project root.
    execute_command(command, PROJECT_ROOT, "build.log", {})
   
 
 
 
# ==============================================================================
# MAIN COMMAND-LINE INTERFACE PARSER
# ==============================================================================
 
def _add_run_parser(subparsers):
    """Attach `run` parser with staged execution and dry-run support."""
    p_run = subparsers.add_parser(
        "run",
        help="Execute a simulation workflow (solve and/or post-process).",
        formatter_class=argparse.RawTextHelpFormatter,
        description=(
            "Execute solver and/or post-processing stages.\n\n"
            "Examples:\n"
            "  pic.flow run --solve --post-process -n 8 --case case.yml --solver solver.yml --monitor monitor.yml --post post.yml\n"
            "  pic.flow run --post-process --run-dir runs/my_run --post post.yml\n"
            "  pic.flow run --solve --post-process --case case.yml --solver solver.yml --monitor monitor.yml --post post.yml --cluster cluster.yml --no-submit\n"
            "  pic.flow run --solve --post-process --case case.yml --solver solver.yml --monitor monitor.yml --post post.yml --dry-run\n"
            "  pic.flow run --solve --post-process --case case.yml --solver solver.yml --monitor monitor.yml --post post.yml --dry-run --format json"
        ),
        epilog="Next: run `pic.flow validate ...` first for config-only checks.",
    )
    run_group = p_run.add_argument_group("stages")
    run_group.add_argument("--solve", action="store_true", help="Execute the solver stage (creates a new run directory).")
    run_group.add_argument("--post-process", action="store_true", help="Execute the post-processing stage on a run directory.")
 
    solver_group = p_run.add_argument_group("solver inputs (required for --solve)")
    solver_group.add_argument("--case", help="Path to the case definition file (e.g., case.yml).")
    solver_group.add_argument("--solver", help="Path to the solver settings profile (e.g., solver.yml).")
    solver_group.add_argument("--monitor", help="Path to the monitoring and I/O profile (e.g., monitor.yml).")
 
    post_group = p_run.add_argument_group("post-processor inputs (required for --post-process)")
    post_group.add_argument("--run-dir", help="Path to an existing run directory to post-process.\n(Not needed if running with --solve in the same command).")
    post_group.add_argument("--post", help="Path to the post-processing recipe file (e.g., post.yml).")
 
    p_run.add_argument("-n", "--num-procs", type=int, default=1, help="Number of MPI processes for either stage.")
    p_run.add_argument("--cluster", help="Path to cluster.yml for Slurm execution mode.")
    p_run.add_argument("--scheduler", help="Explicit scheduler selector (currently 'slurm').")
    p_run.add_argument("--no-submit", action="store_true", help="Generate Slurm scripts/manifests but do not call sbatch.")
    p_run.add_argument("--dry-run", action="store_true", help="Resolve and print planned commands/artifacts without writing files.")
    p_run.add_argument(
        "--format",
        dest="output_format",
        choices=["text", "json"],
        default="text",
        help="Output format for --dry-run (default: text).",
    )
    return p_run
 
 
def _add_sweep_parser(subparsers):
    p_sweep = subparsers.add_parser(
        "sweep",
        help="Launch a Slurm-based parameter sweep/study.",
        formatter_class=argparse.RawTextHelpFormatter,
        description=(
            "Launch matrix studies from study.yml + cluster.yml.\n\n"
            "Examples:\n"
            "  pic.flow sweep --study study.yml --cluster cluster.yml\n"
            "  pic.flow sweep --study study.yml --cluster cluster.yml --no-submit"
        ),
        epilog="Next: inspect studies/<study_id>/results/metrics_table.csv for aggregated metrics.",
    )
    p_sweep.add_argument("--study", required=True, help="Path to study.yml defining parameter matrix and metrics.")
    p_sweep.add_argument("--cluster", required=True, help="Path to cluster.yml defining Slurm resources.")
    p_sweep.add_argument("--no-submit", action="store_true", help="Generate all study artifacts without submitting jobs.")
    return p_sweep
 
 
def _add_validate_parser(subparsers):
    p_validate = subparsers.add_parser(
        "validate",
        help="Validate config files without launching solver/post.",
        formatter_class=argparse.RawTextHelpFormatter,
        description=(
            "Validate one or more config roles. No solver/post execution and no run/study artifact writes.\n\n"
            "Examples:\n"
            "  pic.flow validate --case case.yml --solver solver.yml --monitor monitor.yml\n"
            "  pic.flow validate --post post.yml --cluster cluster.yml\n"
            "  pic.flow validate --study study.yml --cluster cluster.yml --strict"
        ),
        epilog="Next: run `pic.flow run --dry-run ...` to inspect resolved commands/artifacts.",
    )
    p_validate.add_argument("--case", help="Path to case.yml")
    p_validate.add_argument("--solver", help="Path to solver.yml")
    p_validate.add_argument("--monitor", help="Path to monitor.yml")
    p_validate.add_argument("--post", help="Path to post.yml")
    p_validate.add_argument("--cluster", help="Path to cluster.yml")
    p_validate.add_argument("--study", help="Path to study.yml")
    p_validate.add_argument("--strict", action="store_true", help="Enable additional strict checks for selected roles.")
    return p_validate
 
 
def _add_init_parser(subparsers):
    p_init = subparsers.add_parser(
        "init",
        help="Initialize a new case study directory from a template.",
        formatter_class=argparse.RawTextHelpFormatter,
        description=(
            "Create a study directory from examples/<template_name>.\n\n"
            "Examples:\n"
            "  pic.flow init flat_channel --dest my_case\n"
            "  pic.flow init bent_channel --dest my_bent_case --copy-binaries"
        ),
        epilog="Next: run `pic.flow validate --case ... --solver ... --monitor ...` before execution.",
    )
    p_init.add_argument("template_name", help="Name of the case template directory to copy (e.g., 'flat_channel').")
    p_init.add_argument(
        "--dest",
        dest="dest_name",
        help="Optional name for the new directory. Defaults to the template name.\nPath is relative to your current working directory.",
    )
    p_init.add_argument(
        "--copy-binaries",
        action="store_true",
        help="Copy executables into the case directory instead of symlinking.\nThis creates a fully portable, self-contained case study.",
    )
    return p_init
 
 
def _add_build_parser(subparsers):
    p_build = subparsers.add_parser(
        "build",
        help="Build project executables using the Makefile.",
        formatter_class=argparse.RawTextHelpFormatter,
        description=(
            "Calls the project's build.sh script. Any arguments provided after 'build'\n"
            "are passed directly to make.\n\n"
            "Examples:\n"
            "  pic.flow build\n"
            "  pic.flow build clean-project\n"
            "  pic.flow build SYSTEM=cluster\n"
            "  pic.flow build postprocessor"
        ),
        epilog="Next: run `pic.flow --help` or `pic.flow run --help` for execution commands.",
    )
    p_build.add_argument(
        "make_args",
        nargs=argparse.REMAINDER,
        help="Arguments to pass directly to the make command (e.g., 'clean-project').",
    )
    return p_build
 
 
def build_main_parser():
    """Build and return the top-level CLI parser."""
    parser = argparse.ArgumentParser(
        description="pic.flow: A comprehensive conductor for the PIC-Flow CFD simulation platform.",
        formatter_class=argparse.RawTextHelpFormatter,
        epilog=(
            "Examples:\n"
            "  pic.flow validate --case case.yml --solver solver.yml --monitor monitor.yml --post post.yml\n"
            "  pic.flow run --solve --post-process --case case.yml --solver solver.yml --monitor monitor.yml --post post.yml --dry-run\n"
            "  pic.flow run --solve --post-process --case case.yml --solver solver.yml --monitor monitor.yml --post post.yml --cluster cluster.yml --no-submit\n"
            "  pic.flow sweep --study study.yml --cluster cluster.yml\n\n"
            "Next commands:\n"
            "  - First run: pic.flow init ... -> pic.flow validate ... -> pic.flow run ...\n"
            "  - Config debugging: pic.flow validate ...\n"
            "  - Launch planning: pic.flow run ... --dry-run"
        ),
    )
    subparsers = parser.add_subparsers(dest="command", required=True, help="Available commands")
    _add_run_parser(subparsers)
    _add_sweep_parser(subparsers)
    _add_validate_parser(subparsers)
    _add_init_parser(subparsers)
    _add_build_parser(subparsers)
    return parser
 
 
def dispatch_command(args):
    """Validate argument combinations and dispatch to command handlers."""
    if args.command == "run":
        if not args.solve and not args.post_process:
            fail_cli_usage("At least one stage (--solve or --post-process) must be selected.")
        if args.solve and (not args.case or not args.solver or not args.monitor):
            fail_cli_usage("--solve requires --case, --solver, and --monitor.")
        if args.post_process and not args.post:
            fail_cli_usage("--post-process requires --post.")
        if args.scheduler and not args.cluster:
            fail_cli_usage("--scheduler requires --cluster in this version.")
        run_workflow(args)
        return
    if args.command == "sweep":
        sweep_workflow(args)
        return
    if args.command == "validate":
        validate_workflow(args)
        return
    if args.command == "init":
        init_case(args)
        return
    if args.command == "build":
        build_project(args)
        return
    fail_cli_usage(f"Unsupported command '{args.command}'.")
 
 
if __name__ == "__main__":
    main_parser = build_main_parser()
    dispatch_command(main_parser.parse_args())