Python驱动Mumax3微磁学仿真：从脚本散落到可复用工具链

简介mx3tools 是一套面向微磁学仿真研究者的 Python 工具集合专为配合 Mumax3 运行与分析仿真而设计适合已具备一定 Python 基础、需要批量处理仿真输出的科研人员与研究生使用。资源包共 20 个文件以 12 个 py 脚本为核心辅以 ovf 数据样例、txt 说明、yaml 与 cfg 配置、md 文档及 license 等压缩包约 8.55MB结构紧凑、便于直接安装调用。其核心能力在于将 simulation.out 目录下的一组 .ovf 文件批量读入 numpy 数组返回形状为 (n_files, nz, ny, nx, 3) 的数据方便后续对磁化分量做统计、可视化与绘图分析同时提供 ioutil、plotutil、statutil、datautil、simrun 等模块覆盖数据读写、绘图、统计与仿真运行等环节。目前已有 1521 人学习下载可作为微磁仿真数据处理流程中的实用脚本参考帮助读者减少重复的格式解析工作把精力集中在物理结果分析上。1. 从脚本散落到工具链mx3tools 到底解决什么如果你用 Mumax3 做过微磁学仿真大概率经历过这个循环手写.mx3脚本、命令行跑一遍、打开 OVF 文件、用 Python 临时拼几行 matplotlib 画图、再手动改参数重跑。单次实验没问题一旦要扫 50 组 DMI 强度或者 20 组几何尺寸脚本就开始散落一地命名靠test_final_v3这种玄学后缀撑着。mx3tools 这个方向要解决的就是把「跑仿真」和「分析仿真」这两件事从一次性脚本变成可复用的 Python 工具链——用 Python 生成 Mumax3 输入、批量调度、解析输出、统一出图。它适合两类人一是刚开始用 Mumax3、想少走弯路的研究生二是已经跑了几百次仿真、被文件管理折磨过的从业者。下面我按自己搭这套流程的顺序把选型、代码、参数和踩过的坑讲清楚。2. 用 Python 驱动 Mumax3输入生成与批量调度2.1 为什么不让 Python 直接调 Mumax3 二进制Mumax3 本身是个 Go 写的命令行程序输入是.mx3文本脚本输出是table.txt、*.ovf、*.png这些文件。它没有官方 Python API所以「用 Python 运行仿真」的本质是进程调度 文件读写不是库调用。常见做法有两层第一层是模板化输入。把.mx3里会变的参数交换刚度、DMI 系数、外场、几何尺寸抽成占位符Python 用字符串替换或 Jinja2 渲染出具体脚本。第二层是子进程调度。用subprocess调mumax3可执行文件配合concurrent.futures做并发或者干脆生成一批脚本丢给集群的作业系统。我一般会先确认环境里mumax3在 PATH 上且 GPU 可用# 确认 mumax3 可执行文件位置和版本 which mumax3 mumax3 -version # 确认 GPU 状态Mumax3 依赖 CUDA nvidia-smiwhich找不到就说明没配 PATH后面subprocess会直接报FileNotFoundError。nvidia-smi看不到 GPU 或者驱动版本太老Mumax3 启动时会报 CUDA 初始化失败这种错误在批量跑的时候会淹没在日志里所以先单独确认一次。2.2 一个可复用的输入模板与渲染函数下面这个模板覆盖了最常见的层状结构仿真一个矩形薄膜加 DMI 和垂直各向异性扫外场。占位符用{{ }}标出方便替换。# mx3_template.py # 一个参数化的 Mumax3 输入模板用 Python 渲染成具体 .mx3 文件 TEMPLATE setgridsize({{ nx }}, {{ ny }}, 1) setcellsize({{ cx }}e-9, {{ cy }}e-9, {{ cz }}e-9) Msat {{ msat }}e6 Aex {{ aex }}e-12 Dind {{ dmi }}e-3 Ku1 {{ ku1 }}e6 AnisU vector(0, 0, 1) alpha {{ alpha }} B_ext vector({{ bx }}e-3, 0, {{ bz }}e-3) m uniform(1, 0, 0) relax() tableadd(E_total) tableadd(m) autosave(m, {{ save_every }}e-12) run({{ runtime }}e-9) def render_mx3(params: dict) - str: 把参数字典渲染成 Mumax3 脚本字符串。 params 必须包含模板里所有占位符对应的键。 return TEMPLATE.format(**params) if __name__ __main__: p dict(nx128, ny128, cx4, cy4, cz4, msat0.8, aex13, dmi1.5, ku10.4, alpha0.02, bx0, bz50, save_every50, runtime2) print(render_mx3(p))逻辑说明TEMPLATE用 Python 的str.format做替换占位符名和params的键一一对应。render_mx3只做渲染不碰文件系统这样单元测试时可以直接断言输出字符串。参数说明nx/ny是网格点数cx/cy/cz是单胞尺寸纳米msat单位是 1e6 A/maex单位 1e-12 J/mdmi单位 1e-3 J/m²ku1单位 1e6 J/m³bz是外场mT。这些单位前缀是 Mumax3 的惯例写错量级是新手最常见的翻车点——比如把aex13当成 13 J/m 而不是 13e-12仿真结果会离谱到没法看。2.3 批量调度并发跑多个参数点有了渲染函数批量就是「生成脚本 → 建目录 → 起进程」。关键点是每个参数点一个独立工作目录避免输出文件互相覆盖。# batch_run.py import os, subprocess, itertools from concurrent.futures import ThreadPoolExecutor, as_completed from mx3_template import render_mx3 def run_one(workdir: str, params: dict, mumax_binmumax3): os.makedirs(workdir, exist_okTrue) mx3_path os.path.join(workdir, sim.mx3) with open(mx3_path, w) as f: f.write(render_mx3(params)) # -f 覆盖已有输出-o 指定输出目录 cmd [mumax_bin, -f, -o, workdir, mx3_path] r subprocess.run(cmd, capture_outputTrue, textTrue) return workdir, r.returncode, r.stderr[-500:] def sweep(dmi_values, bz_values, max_workers2): jobs [] for dmi, bz in itertools.product(dmi_values, bz_values): tag fdmi{dmi:.2f}_bz{bz:.0f} params dict(nx128, ny128, cx4, cy4, cz4, msat0.8, aex13, dmidmi, ku10.4, alpha0.02, bx0, bzbz, save_every50, runtime2) jobs.append((os.path.join(runs, tag), params)) results [] with ThreadPoolExecutor(max_workersmax_workers) as ex: futs [ex.submit(run_one, wd, p) for wd, p in jobs] for fut in as_completed(futs): results.append(fut.result()) return results逻辑说明run_one负责单点sweep用itertools.product生成参数组合ThreadPoolExecutor控制并发。参数说明max_workers不要设太大Mumax3 吃满 GPU同一张卡上并发两个以上通常只会互相拖慢甚至显存不够直接崩。我一般设 1 到 2具体看单次仿真的显存占用。-f是强制覆盖-o指定输出目录这两个参数不加的话输出会散在当前目录批量跑完根本分不清哪个文件属于哪个参数点。r.stderr[-500:]只留最后 500 字符是因为 Mumax3 的日志很长全存下来内存吃不消但报错信息通常在末尾。3. 解析 OVF 与 table.txt把仿真输出变成数组3.1 table.txt 的读取与列名陷阱table.txt是 Mumax3 的标量时间序列输出格式是带表头的空白分隔文本。看起来简单但列名会随tableadd的调用顺序变化硬编码列索引迟早出事。# parse_table.py import numpy as np def load_table(path: str) - dict: 读取 Mumax3 的 table.txt返回 {列名: np.ndarray}。 自动跳过以 # 开头的注释行用第一行非注释行做表头。 header None rows [] with open(path) as f: for line in f: line line.strip() if not line or line.startswith(#): continue if header is None: header line.split() continue rows.append([float(x) for x in line.split()]) data np.array(rows) return {name: data[:, i] for i, name in enumerate(header)}逻辑说明逐行读遇到#跳过第一个非注释行当表头之后的行转 float。参数说明返回的是字典而不是二维数组这样调用方用t[E_total]而不是t[:, 3]列顺序变了也不影响。注意 Mumax3 的表头里可能有m_x、m_y、m_z这种带下划线的名字用字典键访问最稳。如果某次仿真中途崩了table.txt最后一行可能是半截的float()会抛ValueError所以生产环境里最好包一层 try 跳过坏行。3.2 OVF 矢量场的读取自己写还是用库OVFOpen Visualization Format是 Mumax3 存磁化矢量场的主要格式二进制或文本都有。常见做法是用现成的discretisedfield库它和 Mumax3 同源读 OVF 很顺。但如果你的环境不方便装这个库自己解析文本版 OVF 也不难。# parse_ovf.py import numpy as np def load_ovf_text(path: str): 读取文本格式的 OVF 1.0 文件返回 (data, meta)。 data 形状为 (nx, ny, nz, 3)。仅支持文本段二进制段会报错。 meta {} with open(path, rb) as f: raw f.read() # 头部是 # 开头的键值对直到遇到数据段标记 text raw.decode(latin-1) lines text.splitlines() data_start None for i, line in enumerate(lines): if line.startswith(#): if : in line: k, v line[1:].split(:, 1) meta[k.strip()] v.strip() else: data_start i break if meta.get(Begin, ).strip() ! Segment: raise ValueError(未找到数据段起始标记) # 数据段每行一个值按 x 最快、z 最慢的顺序排列 values [] for line in lines[data_start 1:]: line line.strip() if not line or line.startswith(#): continue values.append(float(line)) arr np.array(values) nx int(meta[xnodes]); ny int(meta[ynodes]); nz int(meta[znodes]) arr arr.reshape(nz, ny, nx, 3) # OVF 是 z 最慢 return np.transpose(arr, (2, 1, 0, 3)), meta逻辑说明先按#解析头部元数据找到数据段后逐行读浮点数再按xnodes/ynodes/znodesreshape。参数说明OVF 的存储顺序是 x 变化最快、z 最慢所以 reshape 成(nz, ny, nx, 3)后再转置成(nx, ny, nz, 3)这样data[i, j, k]对应网格点(i, j, k)的磁化矢量。latin-1解码是为了兼容某些带非 ASCII 注释的文件。如果文件是二进制 OVF这个函数会在float()处报错那就得换discretisedfield或者按二进制格式另写解析。我一般优先用discretisedfield自己写只在依赖受限时用。3.3 从矢量场算标量拓扑荷与磁化分量拿到(nx, ny, nz, 3)的数组后最常用的分析是算磁化分量平均值和拓扑荷skyrmion number。拓扑荷的离散计算要注意边界处理。# analyze.py import numpy as np def magnetization_components(m): 返回 mx, my, mz 三个分量数组。 return m[..., 0], m[..., 1], m[..., 2] def skyrmion_number(m): 用离散格点法计算拓扑荷。 m 形状 (nx, ny, nz, 3)取 nz0 的切片做二维计算。 mx m[:, :, 0, 0]; my m[:, :, 0, 1]; mz m[:, :, 0, 2] # 归一化防止数值误差导致模长偏离 1 norm np.sqrt(mx**2 my**2 mz**2) mx, my, mz mx/norm, my/norm, mz/norm # 相邻格点的叉积和点积累加立体角 def cross_dot(a, b, c): return np.sum(a * np.cross(b, c), axis-1) mvec np.stack([mx, my, mz], axis-1) # 四个角点构成一个 plaquette m00 mvec[:-1, :-1] m10 mvec[1:, :-1] m11 mvec[1:, 1:] m01 mvec[:-1, 1:] omega (cross_dot(m00, m10, m11) cross_dot(m00, m11, m01)) # 用 arctan2 得到每个 plaquette 的立体角 area np.arctan2(omega, 1 np.sum(m00*m10, -1) np.sum(m10*m11, -1) np.sum(m11*m01, -1) np.sum(m01*m00, -1)) return np.sum(area) / (4 * np.pi)逻辑说明拓扑荷本质是磁化场在单位球面上覆盖的立体角除以 4π。这里用四个相邻格点构成 plaquette算叉积点积组合再取arctan2比直接求和更稳定。参数说明m必须是归一化后的单位矢量如果仿真输出的磁化模长有微小偏差先归一化再算否则拓扑荷会漂。取nz0切片是因为二维 skyrmion 通常只在单层里算多层结构要对每层分别算再求和。这个函数返回的是无量纲数理想 skyrmion 接近 ±1数值误差在 1e-3 量级算正常。4. 避坑与排查批量仿真里最容易翻车的五件事4.1 现象批量跑完发现一半目录是空的原因subprocess.run没检查返回码Mumax3 启动失败比如显存不足时进程直接退出但 Python 侧没抛异常目录建了、脚本写了、输出没有。解决run_one里必须判断r.returncode ! 0把 stderr 落盘到error.log批量结束后统一检查哪些点失败。我一般会在sweep返回结果里过滤出失败项重跑前先看错误日志别盲目重试。4.2 现象table.txt 列名对不上画出来的曲线是错的原因模板里tableadd的调用顺序变了或者不同批次的脚本模板不一致导致列索引偏移。解决永远用列名字典访问不要用data[:, 3]这种硬编码。如果列名本身会变在load_table里加断言检查期望的列名是否存在缺了就早报错。4.3 现象OVF 读出来形状不对图像是花的原因reshape 顺序搞反了。OVF 是 x 最快、z 最慢如果按(nx, ny, nz, 3)直接 reshape 而不转置空间坐标就错位了。解决严格按(nz, ny, nx, 3)reshape 再转置并且用一个小测试文件验证——比如均匀磁化m(1,0,0)读出来应该每个点都是(1,0,0)形状对了值不对就是顺序问题。4.4 现象并发跑多个仿真速度反而比串行慢原因Mumax3 吃满 GPU多个进程抢同一张卡上下文切换和显存竞争让总时间变长甚至触发 OOM。解决max_workers设 1靠参数点数量堆吞吐而不是靠并发。如果确实要并发先测单次仿真的显存占用确保max_workers × 单次显存 总显存并且留 20% 余量。多卡环境可以用CUDA_VISIBLE_DEVICES给每个进程绑不同卡。4.5 现象拓扑荷算出来是 0.3 这种非整数原因磁化场没有完全弛豫到稳态或者边界格点的磁化被钉扎导致立体角积分不闭合。解决先确认relax()收敛了看table.txt里能量是否平稳再检查边界——如果边界磁化不是均匀朝上拓扑荷本来就不该是整数。另外网格太粗比如 32×32时离散误差大拓扑荷会偏离加密到 128×128 以上通常能收敛到 ±1 附近。5. 进阶把分析流程封装成可复现的管线前面几章是零件这一章讲怎么拼成一条能反复用的管线。核心思路是配置驱动 结果落盘 幂等重跑。我一般会用一个 YAML 或 JSON 描述整个扫描计划Python 读配置、生成任务、跑完把每个点的关键指标能量、拓扑荷、磁化分量均值写进一个汇总 CSV。这样换台机器、换个人只要配置和代码在结果就能复现。# pipeline.py import json, os, csv from batch_run import sweep from parse_table import load_table from parse_ovf import load_ovf_text from analyze import skyrmion_number def summarize(workdir: str) - dict: 从一个仿真目录里提取关键指标。 out {workdir: workdir} table_path os.path.join(workdir, table.txt) if os.path.exists(table_path): t load_table(table_path) out[E_final] float(t[E_total][-1]) # 找最后一个 OVF 快照 ovfs sorted(f for f in os.listdir(workdir) if f.endswith(.ovf)) if ovfs: m, _ load_ovf_text(os.path.join(workdir, ovfs[-1])) out[sk_number] float(skyrmion_number(m)) return out def run_pipeline(config_path: str, out_csvsummary.csv): with open(config_path) as f: cfg json.load(f) results sweep(cfg[dmi_values], cfg[bz_values], max_workerscfg.get(max_workers, 1)) rows [] for workdir, rc, err in results: if rc ! 0: rows.append({workdir: workdir, error: err}) continue rows.append(summarize(workdir)) # 汇总写 CSV字段取并集 keys set() for r in rows: keys.update(r.keys()) with open(out_csv, w, newline) as f: w csv.DictWriter(f, fieldnamessorted(keys)) w.writeheader() w.writerows(rows) return rows逻辑说明summarize从单个目录抽指标run_pipeline读配置、跑扫描、汇总。参数说明config.json里放dmi_values、bz_values、max_workers这样扫描范围不用改代码。summarize里对table.txt和 OVF 都做了存在性检查因为失败的目录可能只有部分输出。汇总 CSV 用字段并集避免某个点缺字段导致DictWriter报错。一个具体技巧给每个参数点算一个稳定的哈希标签而不是用dmi1.50_bz50这种浮点字符串。浮点格式化在不同机器上可能差一位导致目录名不一致、重跑时对不上。我一般用hashlib.md5(json.dumps(params, sort_keysTrue).encode()).hexdigest()[:8]做标签配置一样标签就一样幂等重跑时直接跳过已完成的目录。验证方法上我习惯先跑一个已知解析解或已知文献值的基准点——比如无 DMI、无外场时应该得到均匀磁化拓扑荷为 0或者一个已知 skyrmion 稳定窗口的参数点拓扑荷接近 1。基准点跑通再铺开扫描否则批量跑完发现全错返工成本很高。最后说个血泪经验别在仿真跑完的当天就删中间文件。我有一次扫完 200 个点觉得 OVF 太占地方就清了结果第二天发现拓扑荷算法有个边界 bug想重算却没有原始场数据只能全部重跑。现在我的习惯是 OVF 至少留一周或者先压缩归档再删。这套流程本身不复杂难的是忍住「先跑起来再说」的冲动把配置、日志、基准点这几步做扎实。希望帮到你。本文还有配套的精品资源点击获取