Z-mag：一种用于有限元微磁模拟的双图立体坐标表示方法

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：Z-mag: A two-chart stereographic formulation for finite-element micromagnetic simulations

　　•Z-mag：基于双图立体投影的有限元微磁学算法；通过构造确保|m|=1。•使用符号参数化的对偶图表并进行连续加权，以实现平滑的演化过程。•采用算子分割的LLG方法、牛顿-克里洛夫（Newton-Krylov）算法、有限元/边界元（FEM/BEM） Magnetostatics方法以及热力学控制。

引言
微磁学模拟在理解和预测纳米尺度材料的磁行为中起着核心作用，其应用范围涵盖自旋电子学、传感器、磁存储设备和量子计算[1]、[2]、[3]、[4]、[5]、[6]。微磁学采用连续介质近似方法，在介观尺度上对离散的原子磁矩进行空间平均，以获得连续的磁化矢量场[1]、[7]。在远低于居里温度的情况下，磁化强度保持饱和状态，磁化动力学仅由重定向过程控制，从而产生满足点态约束|m|=1的单位矢量场m。这种动力学过程由朗道-利夫希茨-吉尔伯特（LLG）方程描述，该方程决定了在外加磁场作用下的进动和阻尼运动。主要的计算挑战在于在离散层次上求解LLG方程的同时保持单位长度约束。

传统的微磁学建模方法主要基于两种参数化方式：笛卡尔坐标系和球坐标系[8]、[9]。笛卡尔坐标系将磁化表示为m=(mx,my,mz)，并通过归一化、惩罚函数或拉格朗日乘数来强制满足单位长度约束[1]、[10]、[11]、[12]、[13]。在基于归一化的方法中，磁化值会定期重新缩放到单位长度，即m←m/|m|。这可以在选定的离散状态下保证|m|=1，但在中间迭代过程中该约束并不被内在地执行，因此长时间模拟后可能会产生漂移。在低阻尼或弱耗散情况下，重复归一化可能会引入人为的耗散和相位误差，从而扰乱物理进动过程。惩罚函数通过添加形式为λ(|m|2?1)2的项来更平滑地执行约束，其中λ>0是惩罚参数。λ的选择很微妙：如果λ太小，则约束效果较弱；如果太大，则离散问题的条件数会变得非常差。拉格朗日乘数方法引入了一个辅助场μ，并通过积分Ωmμ(|m|2?1)dV来弱化约束，其中Ωm?R3表示磁域。虽然这种方法在连续层面上可以精确执行约束，但它会增加未知数的数量，并导致系统变为鞍点问题，通常需要使用inf-sup-stable离散化技术和专门的求解器或预处理器。尽管存在这些缺点，但由于其概念上的简单性和实现的便利性，笛卡尔坐标系仍被广泛使用。

球坐标系通过极角θ0∈[0,π]和方位角θ1∈[0,2π)来参数化磁化，即m=(sinθ0cosθ1,sinθ0sinθ1,cosθ0)，从而通过构造满足单位长度约束[14]、[15]、[16]、[17]、[18]。其主要缺点是在极点θ0=0和θ0=π处存在坐标奇异性，此时θ1的几何意义不明确。在这些点上，度量因子sin2θ0消失，因此θ1的变化不会在单位球面S2上产生物理位移。在离散化表示中，这会导致与方位角方向相关的算子条件数严重恶化或秩不足，如果没有专门的处理方法，可能会导致求解器停滞或无法收敛。实际实现中通常采用避免极点的措施，例如角度限制或基于插值的正则化[18]、[19]。尽管这些措施可以有效，但它们会增加算法复杂性，并且在磁化值接近或穿过极点时可能会引入伪影。这限制了球坐标系在通用模拟中的鲁棒性，特别是在涉及强垂直磁各向异性系统、涡旋核心或穿过极点的磁畴壁的情况下。

立体投影提供了一种替代笛卡尔和球坐标系的方法，它将单位球面S2保形映射到平面R2上。给定一个投影极点P，S2内除了P之外的任何状态m都可以通过直线P和m与投影平面的交点q=(u,v)来表示。通过构造确保单位长度约束|m|=1，而保形度量使得LLG方程的变分形式中的切平面具有各向同性，避免了球坐标系的方位角退化问题。然而，当m→P时，坐标会发散（|q|→∞），导致算子在投影极点附近变得刚性增大，图表分辨率降低。之前的工作已经在特定场景中解决了这一问题：Horley等人[9]、[20]使用了基于阈值的极点切换方法来处理均匀磁化动力学，而Rybakov等人[21]、[22]在手性纹理的能量最小化过程中采用了自适应图表重新定义方法。这些方法依赖于离散的切换事件，并未提供适用于时变微磁学和耦合磁静力学的通用有限元框架。

在数值实现方面，微磁学软件主要基于两种空间离散化方法。有限差分（FD）方法，如OOMMF[23]和MuMax3[24]，利用规则网格和快速傅里叶变换（FFT）来快速计算磁静力场，但通常仅限于直线几何形状和均匀材料参数。有限元（FE）方法，如magnum.fe[25]，提供了更大的几何和材料灵活性，包括曲线边界、异质介质和局部网格细化，但代价是计算和内存需求更高，因为使用了非结构化离散化和密集的退磁算子。在这两种方法中，大多数商用代码都使用带有显式范数控制的笛卡尔磁化分量。据我们所知，主流通用微磁学软件包中并未使用球坐标或立体坐标作为主要参数化方式。

在这项工作中，我们介绍了Z-mag，这是一个开源的有限元微磁学求解器，它基于磁化的双图立体投影表示。该方法演化了两个互补的立体场qA和qB，它们与对偶极点相关，并通过离散算子中的空间变化权重将它们结合起来。与早期研究中的基于切换的立体投影方法[9]、[20]、[21]、[22]不同，Z-mag在每个时间步长内通过重叠加权算子同时推进两个图表，并通过接受后关闭/重建来确保跨图表的一致性。这种重叠加权结构使得在实际应用中保持有效的图表表示，并避免了时间步进过程中的硬性全局图表切换。立体场和磁静力势在混合FE设置中被离散化。每个时间步长通过交替进行磁化和磁静力求解的算子分割工作流来推进：磁化子问题通过牛顿-克里洛夫迭代处理，而磁静力则通过有限元/边界元（FEM/BEM）方法更新。时间积分使用自适应步长向后微分公式，并结合热力学接受标准。Z-mag建立在FEniCSx/DOLFINx[26]和PETSc[27]的基础上，为复杂几何形状和一般磁化状态提供了稳健的计算框架。该软件以开源形式发布[28]。

本文的其余部分组织如下：第2节介绍双图立体坐标公式和LLG方程；第3节描述有限元实现，包括图集加权、混合离散化、算子分割和自适应时间积分；第4节详细说明图集管理、重叠权重和图表同步；第5节通过基准测试验证求解器的性能；第6节总结我们的贡献。

数学公式
本节建立了用于微磁动力学的双图立体坐标公式的数学基础。在回顾了标准的微磁能量泛函和笛卡尔形式的朗道-利夫希茨-吉尔伯特（LLG）方程之后，我们指出了离散时间积分中的约束执行问题，这些问题促使人们寻求替代参数化方式。然后逐步发展了立体坐标框架：首先研究了球坐标系中的极点奇异性问题。

有限元实现
在第2节介绍了统一的双图立体坐标公式之后，我们现在在Z-mag中展示了其有限元实现。计算模型在整个域中将两个图表场qA和qB以及磁静力势φ作为同时未知数进行处理。图集权重决定了组装算子中每个元素的局部图表贡献，而接受后的图表关闭确保了两种坐标表示之间的一致性。

图集管理
第3.1节介绍了控制每个元素中哪个图表参与物理计算的图集权重函数wA(K)和wB(K)。本节描述了这些权重在时间演化过程中的管理方式以及在整个模拟过程中保持图表同步的机制。

基准问题
为了评估求解器及其双图立体坐标公式的准确性和鲁棒性，我们考虑了两个基准问题，这两个问题分别探究了动态过程的互补方面。第一个基准问题的研究是在恒定外场下研究空间均匀磁化的进动和阻尼现象。这提供了一个行为已知的良好控制环境，可以直接验证陀螺进动和阻尼现象，以及磁化通过极点时图集的稳定性。

结论
本文介绍了Z-mag，这是一个基于统一双图立体投影公式的微磁学求解器，适用于动态LLG模拟。这种双重立体投影方法解决了笛卡尔坐标系和球坐标系中的坐标奇异性问题。与需要通过归一化或拉格朗日乘数显式执行约束的笛卡尔方法不同，立体坐标系通过构造本身就满足了单位长度约束。

作者贡献声明
Stella Brach：撰写——原始草稿、可视化、验证、监督、软件开发、方法学研究、概念化。Federico Soldà：撰写——审阅与编辑、软件开发、概念化。Kathrin Aebischer：撰写——审阅与编辑、概念化。Mark A. Lantz：撰写——审阅与编辑、监督。

利益冲突声明
作者声明没有已知的财务利益关系或个人关系可能影响本文的报告内容。