Unable to solve problem via NonlinearVariationalProblem (Runtime Error)

variational formulation
1

Hi everyone,

I am currently trying to solve a Nonlinear problem involving the heating of a 2D block in vacuum via a Guassian heat source with no Dirichlet boundary conditions. The aim is to create a steady state temperature map of the block. When I run the NonlinearVariationalSolver, I gain this error message:

Traceback (most recent call last):
  File "test2D.py", line 121, in <module>
    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 nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

I understand that this is a runtime problem as the problem did not converge, I am not sure if the weak form is just not solvable or whether I should use a different solver.
My code is shown below:

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
import math

time = 100.0            # final time
num_steps = 100     # number of time steps
dt = time / num_steps # time step size


#build rectangular mesh to model 2D slice
nx = 32
ny = 16
x0 = 0
y0 = 0
x1 = 4
y1 = 2
pi = 3.141592653589793
T_am = Constant("300")
T_a4 = T_am**4
epsilon = 1  # material emissivity
sigma = 5.67E-8 # W/(m2.K4)
es = epsilon*sigma
Q = Expression('2*A*t*Pmax/(pi*w*w)*exp(-2*pow(x[0]-2, 2)/(w*w))',degree=2, A=1, Pmax=200,w=1,x0=0, y0=0,t=0)
#----------THERMAL--PROPERTIES----------
kappa = 94         #conductivity [W/m K] 
c = 0.444          #Specific Heat Capacity [J/gK]
rho = 7            #Density [g/cm^3]
const = kappa /(c * rho)
#---------------------------------------
 
sw = Point (x0, y0)
ne = Point (x1,y1)

mesh = RectangleMesh(sw, ne, nx, ny, diagonal= "right")
V = FunctionSpace(mesh, 'P', 2)

# Define boundary condition
# build expression for power guassian (x0 and y0 displace the beam from the centre, A is the absorptivity of the target, Pmax is the maximum power of the beam, w is the width of the beam.)


#-----BOUNDARY-CONDITIONS--------------------------


#-----Define Markers for the different parts of the boundary-------
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
tol = 1E-14
class BoundaryX0(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], x0, tol)
    
bx0 = BoundaryX0()
bx0.mark(boundary_markers, 0)

class BoundaryX1(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], x1, tol)
    
bx1 = BoundaryX1()
bx1.mark(boundary_markers, 1)

class BoundaryY0(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], y0, tol)
    
by0 = BoundaryY0()
by0.mark(boundary_markers, 2)

class BoundaryY1(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], y1, tol)
    
by1 = BoundaryY1()
by1.mark(boundary_markers, 3)


# For the implimentation of these boundary conditions to be general, we can let the user specify what kind of boundary condition that applies to each of the four boundaries. We set up a Python dictionary for this purpose.


#ENCODE BOUNDARY CONDITIONS (MAY NEED ALTERING OR REMOVAL)
boundary_conditions = {0: {'Robin':     (es, T_a4)},
                       1: {'Robin':     (es, T_a4)},
                       2: {'Robin':     (es, T_a4)},
                       3: {'Robin':     (es, T_a4)}}


ds = Measure('ds', domain=mesh, subdomain_data = boundary_markers)
integrals_R = []
T = TrialFunction(V)
v = TestFunction(V)
T_ = Function(V)
T_n = interpolate(T_am, V)

#Account for Radiation
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        integrals_R.append(es*(T**4 - T_a4)*v*dt*ds(i))

#----------------------------------------
vtkfile = File('test2D/solution.pvd')
# fix F
F = T*v*dx - const*dt*dot(grad(T), grad(v))*dx + const*Q*v*dt*ds(3) - const*sum(integrals_R) -T_n*v*dx
T_= Function(V)   # the most recently computed solution
F = action(F,T_)
J = derivative(F, T_,T)
t = dt
bcs = []
while t <= time:
    print ('time =', t)
    Q.t = t
    problem = NonlinearVariationalProblem(F, T_, bcs, J)
    solver = NonlinearVariationalSolver(problem)
    solver.solve()
    prm = solver.parameters
    prm['newton_solver']['absolute_tolerance'] = 1E-8
    prm['newton_solver']['relative_tolerance'] = 1E-7
    prm['newton_solver']['maximum_iterations'] = 100
    prm['newton_solver']['relaxation_parameter'] = 1.0
    t += dt
    T_n.assign(T_)
    vtkfile << (T_, t)
    plot(T_)

p = plot(T_)
p.set_cmap("inferno")
#p.set_clim(0,100)
plt.colorbar(p)
plt.savefig("test2D.pdf")

Any help would be appreciated. Thanks again.

2

The error message tells you this:

To start you might add set_log_level(20) at the beginning of your code in order to monitor the Newton iterations. Maybe it is even converging but 100 iterations is not enough?

Smaller time steps could help, too, and ‘continuation’: start with very weak forcing where the solver converges, then increase gradually, using the previous solution as initial guess for Newton, until you reach the quantities/forcings of the original problem.
In any case, I’d try to reduce the complexity/nonlinearity of the problem and find a simpler setting where the solver converges. Then get more complex from there.

Cheers