# NonlinearVariationalSolver converges with zero iteration

417
views
0
7 months ago by
In my quest to implement the two-temperature model (TTM), I have revised my math and (hopefully) progressed with the code compared to my last question (see https://www.allanswered.com/post/jrplb/how-to-implement-the-two-temperature-model/).

I have used NonlinearVariationalSolver as the solver, but what I get is the output "Out[9]: (0, 1)" which, I gather, means it successfully converged with 0 iteration. What can I possibly do to debug and improve this issue?

EDIT: I would like to provide some more details about my problem:

I am using FEniCS as a Docker image. The image id is 822b66610353.

I have previously implemented in FEniCS the heat equation with the same source term f, subdomains and boundaries. I have reused that code to implement the TTM. Therefore, I believe these elements of the code are sound.

One difficulty I had was to come up with the weak form. It is likely to be one of the source of error. I have reviewed it, and this is what I get nonetheless.

Things I have tried, mostly random shots in the dark, include:
- Increasing the number of cells in the mesh ;
- Changing the type of FiniteElement ;
- Adjusting the tolerance tol ;
- Increasing the number of iterations of the nonlinear solver ;
- Adjusting the relative tolerance of the nonlinear solver ;

Here’s the log output when running the code with set_log_level = DEBUG:
Solving nonlinear variational problem.
 Elapsed wall, usr, sys time: 0.000175, 0, 0 (Build sparsity)
 Elapsed wall, usr, sys time: 0.0001618, 0, 0 (Apply (PETScVector))
 Elapsed wall, usr, sys time: 0.0020336, 0, 0 (Init tensor)
 Elapsed wall, usr, sys time: 2.4e-05, 0, 0 (Delete sparsity)
 Elapsed wall, usr, sys time: 3.53e-05, 0, 0 (Apply (PETScVector))
 System assembler: boundary condition 0 applies to axis 0
 Elapsed wall, usr, sys time: 3.33e-05, 0, 0 (DirichletBC init facets)
 Elapsed wall, usr, sys time: 0.0026449, 0, 0 (DirichletBC compute bc)
 Elapsed wall, usr, sys time: 2.31e-05, 0, 0 (Apply (PETScVector))
 Elapsed wall, usr, sys time: 0.018517, 0.01, 0 (Assemble system)
 Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
 Newton solver finished in 0 iterations and 0 linear solver iterations.

Any help or pointers would be appreciated, thanks in advance for your time.
----

Finally, here is the code:
from fenics import *
import numpy as np

# Coefficients
t_i = 0.               # initial time
t_f = 100e-12            # final time
time_zero = 1e-12
num_steps = 100      # number of time steps
dt = (t_f - t_i) / num_steps # time step size
initial_temperature = 300. # K

C_metal = .129      # heat capacity, gold / J g-1 K-1
# K_metal = 3.18      # thermal conductivity, gold / W cm-1 K-1
rho_metal = 19.3    # density, gold / g cm-3

# Smith & Norris (2001) APL 78:1240
Ce_ = 70e-6            # J cm-3 K-2
keq = 317e-6           # W cm-3 K-1
G   = 2.1e10           # W cm-3 K-1
Cl  = C_metal * rho_metal

# Our experiment
en = 1.27e7          # peak power density of the pulse / W cm-2
tau = .06e-12            # pulse width / s
R = .37823           # reflectivity coefficient
alpha = 8.4197e5     # absorption coefficient / cm-1

# Create mesh and define function space
nx = 16
ny = 16
width = 1e-3
height = 1e-3
mesh = RectangleMesh(Point(0., 0.), Point(width, height), nx, ny)

P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

v = TestFunction(V)

u_D = Constant((initial_temperature, initial_temperature))
u = TrialFunction(V)
u_n = interpolate(u_D, V)

Te, Tl = split(u)
Te_n, Tl_n = split(u_n)
ve, vl = split(v)

# Define boundary condition
tol = 1E-12
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[0], 0, tol) or near(x[0], width, tol) or near(x[1], 0, tol) or near(x[1], height, tol))

boundary = Boundary()
boundary_markers = MeshFunction('size_t', mesh, 1)
boundary_markers.set_all(9999)
boundary.mark(boundary_markers, 0)

bcs = []
bcs.append(DirichletBC(V, u_D, boundary_markers, 0))

# We define two subdomains (metal == 0 and water == 1)
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= height / 2. + tol

class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= height / 2. - tol

# materials = CellFunction('size_t', mesh)
materials = MeshFunction('size_t', mesh, 2)    # NB: 'size_t' is the equivalent to 'uint' in newer Dolfin versions
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)

f = Expression('en * (1-refl) * alpha / heat_cap / rho * exp(-pow((x[0]-bw)/w, 2)) * exp(-pow(((t-t0)/tau), 2)) * exp(-alpha*(-x[1]+bh))',
degree=2, w=0.1, tau=tau/2., alpha=alpha, t=0., t0=time_zero, bw=width/2., bh=height/2.,
en=en, refl=R, heat_cap=C_metal, rho=rho_metal)

# Define new measures
dx = Measure('dx', domain=mesh, subdomain_data=materials)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

# Write down the equations
F = Ce_ * Te_n * Te_n * ve * dx \
+ dt * f * ve * dx \
- Ce_ * Te * Te * ve * dx \
- dt * keq * Te / Tl * inner(grad(Te), grad(ve)) * dx \
- dt * G * Te * ve * dx \
+ dt * G * Tl * ve * dx \
+ Cl * Tl * vl * dx \
- dt * G * (Te - Tl) * vl * dx \
- Cl * Tl_n * vl * dx

F_ = action(F, u_n)
J = derivative(F_, u_n, u)

problem = NonlinearVariationalProblem(F_, u_n, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
Best,
François

4
7 months ago by
u_n = interpolate(u_D, V)

bcs.append(DirichletBC(V, u_D, boundary_markers, 0))

The solver has no problem to solve. Everything is in balance.

Also you have to reformulate your weak formulation with dx(0) and dx(1), because as it is, dx get the default value, that is 0. That means that the problem is not solved for the sub-domain marked with 1. The same with ds. Solve the two sub-domains separately with dx(0) and dx(1) and the appropriate boundary conditions on the interface (unless you have regularized constitutive equations).

Also for a mixed problem the boundary conditions are for the subspaces(usually) and not for the mixed-space. So if you want to apply BCs (to the first subspace for example), you replace V with V.sub(0) into the DirichletBC and so on.

Well observed! Indeed, if I set it to something else out of equilibrium, such as u_n = interpolate(Constant((2000, 300)), V), the solver actually iterates once and smooths out u_n to values very close to 300.
written 7 months ago by François Lapointe
I will mark the question as answered, even if your comment now raises issues about the boundary conditions, subdomains and subspaces, issues that I haven't resolved yet. Many thanks!
written 7 months ago by François Lapointe
You should raise another post/question.
written 7 months ago by Minas Mattheakis