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