FiPy模拟铁电双晶晶界特性：200+时间步后图像异常的问题求助

我正在用FiPy模拟铁电双晶的晶界特性，模型基于相场方法，用到了Allen-Cahn和Cahn-Hilliard方程，需要求解四个核心变量：结晶度（eta）、x/y方向的极化强度（对应代码里的ux/uy）以及无量纲电势（volt_hat）。

目前遇到的棘手问题是：用扫掠法结合线性求解器求解偏微分方程时，前200个时间步的结果都完全符合预期，但超过这个步数后，生成的图像就出现了完全不合理的条纹和奇怪图案。我需要让代码能稳定运行更多时间步，以此对标一篇论文的结果，作为我的基准验证步骤，希望能得到大家的帮助。

以下是我的代码片段：

from fipy.solvers import LinearPCGSolver, LinearGMRESSolver, LinearLUSolver
from fipy import MatplotlibViewer
import os
from fipy import dump

# Ensure the directory exists
output_dir = "PickleJar"
os.makedirs(output_dir, exist_ok=True)

# Initialize the solver with a specified tolerance
solver1 = LinearGMRESSolver(tolerance=1e-15)
solver2 = LinearLUSolver(tolerance=1e-10)
print("Solver initialized with tolerance:", solver1.tolerance)

# Other variables
div_u = CellVariable(name='div_u', mesh=mesh)
psi = CellVariable(name="psi", mesh=mesh)
mag_u = CellVariable(name='dimensionless pol', mesh=mesh)

# Save initial state
dump.write(data={
    'crystalline order': eta,
    'dimensionless voltage': volt_hat,
    'theta': theta,
    'globalux': ux,
    'globaluy': uy,
}, filename='PickleJar/IN_phasefullbc14march_Cellvariables_00000.var')

dt = 2 * 5e-5
timesteps = 2500
# timestep_switch = 500
j = 0  # Counter for file naming
save_interval = 500  # Save data every 1000 timesteps

# Initialize viewers for eta, volt_hat, ux, and uy
eta_viewer = MatplotlibViewer(vars=eta, title="eta Field")
volt_hat_viewer = MatplotlibViewer(vars=volt_hat, title="Dimensionless Voltage (volt_hat)")
mag_u_viewer = MatplotlibViewer(vars=mag_u, title="polarization magnitude (mag_u)")
div_u_viewer = MatplotlibViewer(vars=div_u, title="divergence of polarization (div_u)")
psi_viewer = MatplotlibViewer(vars=psi, title="arctan polarization (psi)")

for i in range(timesteps):
    # if i == timestep_switch:
    #     dt = 2 * 5e-4  # increase timestep after 1000 iterations
    #     print(f"\nTimestep increased to {dt} at iteration {i}")
    
    div_u.setValue(ux.grad[0] + uy.grad[1])
    psi.setValue(np.arctan2(uy, ux))
    mag_u.setValue(np.sqrt(ux**2 + uy**2))
    
    # Update old values only once per timestep
    eta.updateOld()
    ux.updateOld()
    uy.updateOld()
    volt_hat.updateOld()
    
    # Plot the current values for eta, volt_hat, ux, and uy
    eta_viewer.plot()
    volt_hat_viewer.plot()
    mag_u_viewer.plot()
    div_u_viewer.plot()
    psi_viewer.plot()

    # Print which timestep we are on
    print(f"\nStarting timestep {i + 1}/{timesteps}")

    # Initialize residual and sweep count
    res1 = res2 = 1e+100
    max_res = 1e+100
    sweep = 0

    # Perform sweeps to achieve convergence within this timestep
    while (max_res := max(res1, res2)) > 1e-3 and sweep < 30:
        # Solve the equations for this timestep
        res1 = etaEq.sweep(dt=dt, solver=solver1)
        res2 = eqn1.sweep(dt=dt, solver=solver1)
        
        sweep += 1

        # Print sweep number and residual
        print(f"  Sweep {sweep}: Residual = {max_res:.5e}")

    # Check if convergence was achieved within the sweep limit
    if max_res <= 1e-3:
        print(f"  Converged after {sweep} sweeps with final residual {max_res:.5e}")
    else:
        print(f"  Max sweeps reached without full convergence. Final residual: {max_res:.5e}")

备注：内容来源于stack exchange，提问作者Kritika Khanal