I'm modeling the magnetization dynamics of a ferrofluid using a custom PDE in COMSOL Multiphysics. The goal is to couple this with the Magnetic Fields interface, using the Langevin model to update the magnetization M over time. My governing PDE for M is of the form: dM/dt = (1/tau) * (M_L - M) where M_L = Ms * (coth(alpha) - 1/alpha) * (H_loc / |H|), and alpha = (mu * |H|) / (kB * T).

I define H_loc as (B / mu0 - M), based on COMSOL's Magnetic Fields interface. The magnetic flux density B is computed from the magnetic vector potential A. I set the constitutive relation as B = mu0 * (H + M) in the "Magnetic Fields" interface.

I'm using a General Form PDE interface to evolve M (with components Mx, My, Mz) and defined source terms as (M_L - M) / tau with proper initial conditions and regularization of |H| to avoid division by zero.

However, the solver keeps failing with errors like:

• "Matrix has zero on diagonal"
• "NaN or Inf found when solving linear system"
• "Failed to find consistent initial values"
• "Last time step is not converged"

My time range is defined as range(0, 1e-9, 1e-6). I've tried tuning mass/damping coefficients and regularizing small terms.

Has anyone successfully implemented a time-dependent nonlinear magnetization model in COMSOL that is tightly coupled with the Magnetic Fields interface? I would appreciate any advice on improving numerical stability or setup strategies.