Problem with residual

Hello, I have the following non-linear problem:

from fenics import *
from mshr import *
from math import *
from dolfin import *


def k(T):
    return 1523.7*(T**(-1.226))

def alpha(lam,T):
    return (4*pi*k(T))/lam
def R(T):
    return (7.188+k(T)**2)/(21.912+k(T)**2)
def C_v(T):
    return 0.6949*exp(2.375*(1e-4)*T)

#parameters

lam = 8e-05
w_spot = 1e-06
l_spot = 2e-06
power = 40
rho = 2.32
A = (2*power)/(pi*w_spot*l_spot)
scan = 10


# Create geometry
b1 = Box(Point(-0.25, -0.5, -0.1), Point(0.25, 0.5, 0))
geometry = b1

# Create mesh
mesh = generate_mesh(geometry, 50)
#mesh = refine(mesh)

# Save to file and plot
File("mesh.pvd") << mesh

V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
def boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, Constant(300.0), boundary)



#Define variational problem
T = Function(V)
v = TestFunction(V)




x = SpatialCoordinate(mesh)
Src = A*(1-R(T))*alpha(lam,T)*exp(alpha(lam,T)*x[2])*exp(-2*(x[0]/w_spot)**2)*exp(-2*(x[1]/l_spot)**2)




F = k(T)*inner(grad(T),grad(v))*dx + Src*v*dx +rho*C_v(T)*scan*T.dx(0)*v.dx(0)*dx

J = derivative(F, T)
problem = NonlinearVariationalProblem(F, T, bc, J)
solver  = NonlinearVariationalSolver(problem)

prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-7
prm['newton_solver']['maximum_iterations'] = 5
prm['newton_solver']['relaxation_parameter'] = 1.0
solver.solve()

The solver does not converge with:
Newton iteration 0: r (abs) = -nan (tol = 1.000e-08) r (rel) = -nan (tol = 1.000e-07)

what could possibly be the reason for this behaviour?

Thanks a lot

As your initial guess for T is zero, k(T) is infinity. Try with a non-zero initial guess.

2 Likes

I tried by interpolating a constant T, but it still gives me the same error

Following @dokken 's advice and changing your line

T = Function(V)

to

T = project(Constant(1), V)

produces the following output:

Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 8.617e+05 (tol = 1.000e-08) r (rel) = 1.000e+00 (tol = 1.000e-07)
Newton iteration 1: r (abs) = 2.009e-10 (tol = 1.000e-08) r (rel) = 2.331e-16 (tol = 1.000e-07)
Newton solver finished in 1 iterations and 1 linear solver iterations.

What did you do exactly?

2 Likes