pcb-afm-simulation-system/torqrippswap-master/torqrippswap/scan_central_mag_arc.py

"""Scan Central Mag Arc [ED] and collect Motor-CAD torque ripple results.

Each run reloads the baseline model before changing only
MagnetCentralArc_HalbachRing. Solver discretization is separately configurable.
"""

from __future__ import annotations

import argparse
import csv
import json
import math
import statistics
import time
import traceback
from pathlib import Path


ROOT = Path(__file__).resolve().parent
DEFAULT_MODEL = ROOT / "MARS-9S8P_SSSR_Halbach_PCB-V1.0.mot"
PARAMETER = "MagnetCentralArc_HalbachRing"
SECTION_PRIORITY = ("E-Magnetics", "Drive", "Losses", "Materials", "Miscellaneous")
METRICS = (
    "Average torque (virtual work)",
    "Torque Ripple (VW)",
    "Torque Ripple (VW) [%]",
    "Cogging Torque Ripple (Vw)",
    "No load speed",
)


def values_inclusive(start: float, stop: float, step: float) -> list[float]:
    if step <= 0 or stop < start:
        raise ValueError("step must be positive and stop must be >= start")
    count = int(math.floor((stop - start) / step + 1e-9))
    values = [round(start + i * step, 10) for i in range(count + 1)]
    if not math.isclose(values[-1], stop, abs_tol=1e-9):
        values.append(float(stop))
    return values


def norm(name: str) -> str:
    return name.split("_(")[0].strip().replace(" ", "").lower()


def parse_export(path: Path) -> dict[str, dict[str, float]]:
    text = None
    for encoding in ("utf-8-sig", "cp1252", "gbk", "latin-1"):
        try:
            text = path.read_text(encoding=encoding)
            break
        except UnicodeDecodeError:
            pass
    if text is None:
        return {}
    result: dict[str, dict[str, float]] = {}
    section = "(root)"
    for raw in text.splitlines():
        line = raw.strip()
        if not line:
            continue
        if ";" not in line:
            section = line
            result.setdefault(section, {})
            continue
        parts = line.split(";")
        try:
            result.setdefault(section, {})[norm(parts[0].strip().strip('"'))] = float(parts[1])
        except (IndexError, ValueError):
            continue
    return result


def pick(result: dict[str, dict[str, float]], key: str):
    wanted = norm(key)
    for section in SECTION_PRIORITY:
        if wanted in result.get(section, {}):
            return result[section][wanted]
    for section in result.values():
        if wanted in section:
            return section[wanted]
    return ""


def existing_ok_values(csv_path: Path) -> set[float]:
    if not csv_path.exists():
        return set()
    with csv_path.open("r", encoding="utf-8-sig", newline="") as fh:
        return {
            float(row["central_mag_arc_ed"])
            for row in csv.DictReader(fh)
            if row.get("status") == "OK"
        }


def write_analysis(csv_path: Path, output: Path) -> None:
    with csv_path.open("r", encoding="utf-8-sig", newline="") as fh:
        rows = [r for r in csv.DictReader(fh) if r.get("status") == "OK"]
    rows.sort(key=lambda r: float(r["central_mag_arc_ed"]))
    valid = [r for r in rows if r.get("torque_ripple_pct")]
    if not valid:
        output.write_text("# Scan analysis\n\nNo successful results to analyze.\n", encoding="utf-8")
        return
    best = min(valid, key=lambda r: float(r["torque_ripple_pct"]))
    baseline = min(valid, key=lambda r: abs(float(r["central_mag_arc_ed"]) - 120.0))
    ripples = [float(r["torque_ripple_pct"]) for r in valid]
    trend = "decreased" if ripples[-1] < ripples[0] else "increased"
    table = [
        "| Central Mag Arc [ED] | Torque Ripple (VW) [%] | Tavg (VW) [Nm] | Status |",
        "|---:|---:|---:|:---|",
    ]
    table += [
        f"| {float(r['central_mag_arc_ed']):g} | {float(r['torque_ripple_pct']):.5g} | "
        f"{float(r['average_torque_vw_nm']):.6g} | {r['status']} |"
        for r in valid
    ]
    text = f"""# Central Mag Arc scan results and analysis

- Successful points: {len(valid)}
- Minimum Torque Ripple: **{float(best['torque_ripple_pct']):.5g}%** at **{float(best['central_mag_arc_ed']):g} EDeg**
- Average torque at the best point: {float(best['average_torque_vw_nm']):.6g} Nm
- Ripple change versus the point nearest 120 EDeg: {float(best['torque_ripple_pct']) - float(baseline['torque_ripple_pct']):+.4g} percentage points
- Across the scanned interval, Torque Ripple {trend}; population standard deviation is {statistics.pstdev(ripples):.4g} percentage points.

Note: Use this scan for relative comparison. Recheck final candidates with 180
points per electrical cycle and a finer airgap mesh.

## Results

{chr(10).join(table)}
"""
    output.write_text(text, encoding="utf-8")


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--model", type=Path, default=DEFAULT_MODEL)
    parser.add_argument("--start", type=float, default=80.0)
    parser.add_argument("--stop", type=float, default=120.0)
    parser.add_argument("--step", type=float, default=1.0)
    parser.add_argument("--torque-points", type=int, default=120,
                        help="Points per electrical cycle; 0 preserves the model value")
    parser.add_argument("--airgap-mesh", type=int, default=960,
                        help="Airgap mesh/layers value; 0 preserves the model value")
    parser.add_argument("--output-dir", type=Path, default=ROOT / "results")
    parser.add_argument("--port", type=int, default=-1,
                        help="Motor-CAD RPC port for connecting to a specific instance")
    parser.add_argument("--new-instance", action="store_true",
                        help="Open a separate Motor-CAD instance")
    parser.add_argument("--no-resume", action="store_true", help="Do not skip successful CSV rows")
    parser.add_argument("--dry-run", action="store_true")
    args = parser.parse_args()

    model = args.model.resolve()
    if not model.exists():
        raise SystemExit(f"Model does not exist: {model}")
    values = values_inclusive(args.start, args.stop, args.step)
    out_dir = args.output_dir.resolve()
    csv_path = out_dir / "central_mag_arc_scan.csv"
    raw_dir = out_dir / "raw"
    analysis_path = out_dir / "central_mag_arc_analysis.md"
    print(f"Model: {model}", flush=True)
    print(f"Scan: {values[0]:g}..{values[-1]:g} EDeg, {len(values)} points; "
          f"TorquePointsPerCycle={args.torque_points or 'model default'}, "
          f"AirgapMesh={args.airgap_mesh or 'model default'}", flush=True)
    print(f"Results: {csv_path}", flush=True)
    if args.dry_run:
        return 0

    out_dir.mkdir(parents=True, exist_ok=True)
    raw_dir.mkdir(parents=True, exist_ok=True)
    completed = set() if args.no_resume else existing_ok_values(csv_path)
    fields = ["central_mag_arc_ed", "torque_ripple_pct", "torque_ripple_nm",
              "average_torque_vw_nm", "cogging_ripple_vw_nm", "no_load_speed_rpm",
              "seconds", "status", "applied_value", "error"]
    new_file = not csv_path.exists() or args.no_resume
    mode = "w" if args.no_resume else "a"

    import ansys.motorcad.core as pymotorcad
    mc = pymotorcad.MotorCAD(port=args.port, open_new_instance=args.new_instance)
    mc.set_visible(True)
    mc.set_variable("MessageDisplayState", 2)
    mc.display_screen("Scripting")
    with csv_path.open(mode, encoding="utf-8-sig", newline="") as fh:
        writer = csv.DictWriter(fh, fieldnames=fields)
        if new_file:
            writer.writeheader()
        try:
            for index, value in enumerate(values, 1):
                if value in completed:
                    print(f"[{index}/{len(values)}] {value:g} EDeg already complete; skipped", flush=True)
                    continue
                started = time.time()
                row = {k: "" for k in fields}
                row.update(central_mag_arc_ed=value, status="FAILED")
                try:
                    mc.load_from_file(str(model))
                    mc.display_screen("Scripting")
                    if args.torque_points:
                        mc.set_variable("TorquePointsPerCycle", args.torque_points)
                    if args.airgap_mesh:
                        mc.set_variable("AirgapMeshPoints_mesh", args.airgap_mesh)
                        mc.set_variable("AirgapMeshPoints_layers", args.airgap_mesh)
                    mc.set_variable(PARAMETER, value)
                    applied = float(mc.get_variable(PARAMETER))
                    row["applied_value"] = applied
                    if not math.isclose(applied, value, abs_tol=1e-6):
                        raise RuntimeError(f"Parameter verification failed: wrote {value}, read {applied}")
                    mc.do_magnetic_calculation()
                    raw_path = raw_dir / f"central_mag_arc_{value:g}.csv"
                    mc.export_results("EMagnetic", str(raw_path))
                    result = parse_export(raw_path)
                    row.update(
                        torque_ripple_pct=pick(result, "Torque Ripple (VW) [%]"),
                        torque_ripple_nm=pick(result, "Torque Ripple (VW)"),
                        average_torque_vw_nm=pick(result, "Average torque (virtual work)"),
                        cogging_ripple_vw_nm=pick(result, "Cogging Torque Ripple (Vw)"),
                        no_load_speed_rpm=pick(result, "No load speed"),
                        status="OK",
                    )
                    if row["torque_ripple_pct"] == "":
                        raise RuntimeError("Torque Ripple (VW) [%] was not found in exported results")
                except Exception as exc:  # noqa: BLE001
                    row["status"] = "FAILED"
                    row["error"] = f"{type(exc).__name__}: {exc}"
                    traceback.print_exc()
                row["seconds"] = round(time.time() - started, 1)
                writer.writerow(row)
                fh.flush()
                print(f"[{index}/{len(values)}] {value:g} EDeg {row['status']}  Ripple={row['torque_ripple_pct']}%  Tavg={row['average_torque_vw_nm']} Nm  {row['seconds']}s", flush=True)
        finally:
            mc.load_from_file(str(model))

    write_analysis(csv_path, analysis_path)
    print(f"Complete: {csv_path}\nAnalysis: {analysis_path}", flush=True)
    return 0


if __name__ == "__main__":
    raise SystemExit(main())