Linear solver can't converge using PETSc Krylov solver

Hello all

Trying to solve a Poisson equation, which depends on the magnetization

\nabla\cdot\nabla\psi=\nabla\cdot\mathbf{M}

where \mathbf{M}=f(\tau)=f(g(\phi)), i.e., the magnetization is dependent on a phase field \phi (using Allen-Cahn for that, and that works OK). Note that the magnetization sharply goes to zero as \phi\rightarrow1. This creates a discontinuous/sharp interface between the magnetic and non-magnetic materials.

I tried all the solvers in dolfin.list_linear_solver_methods() to no success. I’m using version

>>> dolfin.__version__
'2019.1.0'

Any advise? Thanks in advanced!

Here’s the MWE

from dolfin import *
from ufl import tanh

def boundary(x, on_boundary):
  return on_boundary

Lx, Ly = 2.5e-01, 2.5e-01
Nx, Ny =     200,     200

xL, xH = -Lx/2, Lx/2
yL, yH = -Ly/2, Ly/2

mesh = RectangleMesh.create([Point(xL,xH),Point(yL,yH)],[Nx,Ny],\
                             CellType.Type.quadrilateral)
    
V = FunctionSpace(mesh, "Q", 1)

expSt = "(1+tanh(beta*(rad-sqrt(pow(x[0],2)+pow(x[1],2))))/tanh(beta))/2"
  
ph0 = Expression(expSt, beta=20, rad=0.2, degree=1)
phi = interpolate(ph0, V)

T  = 1023.          #temperature
Tc = [1043., 201/3] #Curie temperature
th = 45.            #field angle 

#mag model parameters
b = [0.1098, -0.368, -0.129]

#reduced magnetization
m = lambda t, *b:\
      pow(1-t+DOLFIN_EPS,-b[1])/(1+b[1]*t+b[0]*pow(t,3/2.)+b[2]*pow(t,7/2.))

#reduced temperature
t = T/(Tc[0]*phi + Tc[1]*(1-phi))
M = conditional(gt(t,1),0,m(t,*b))

M = cos(th*pi/180.)*M*Constant((1,0)) + sin(th*pi/180.)*M*Constant((0,1))

#Poisson model
psi  = TrialFunction(V)
psi_ = TestFunction(V)

a = inner(grad(psi),grad(psi_))*dx
L = div(M)*psi_*dx

psi  = Function(V)

sparams = {'linear_solver':'gmres'}
solve(a == L, psi, [], solver_parameters=sparams)

And the error

Solving linear variational problem.
Traceback (most recent call last):
  File "mwe.py", line 52, in <module>
    solve(a == L, psi, [], solver_parameters=sparams)
  File "~/anaconda/3.7/envs/FEniCS3.7/lib/python3.7/site-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "~/anaconda/3.7/envs/FEniCS3.7/lib/python3.7/site-packages/dolfin/fem/solving.py", line 247, in _solve_varproblem
    solver.solve()  
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** 
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
*** 
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset:  f3de96341c54e434c48115de5a3975b482aee933
*** -------------------------------------------------------------------------

The problem is with your mesh generation

which should be:

mesh = RectangleMesh.create([Point(xL, yL), Point(xH, yH)], [Nx, Ny],
                            CellType.Type.quadrilateral)
1 Like

Yes, bad copy/paste. My original issue was with BC definitions. Thanks!