I am trying to simulate grain boundary properties in a ferroelectric bicrystal using FiPy. I have basically three variables to solve namely, degree of crystalinity (eta), polarization magnitude in the x and y direction(polx and poly) and the electric potential (voltage). The problem is based on the phase-field modelling and uses Allen-Cahn and Cahn-Hilliard equations.

While trying to solve the pde's through sweep method, using LinearLUSolver I am facing problem because the code seems to work fine untill 200 time steps, but after this the code doesnot work as expected and the plots are very weird with stripes and patterns I wouldn't imagine. please, help me what I should be doing to make the code work for more number of timesteps because I am trying to simulate the result from a reserach paper which would be my benchmarking step.

The code is given below.

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}")