### A two subdomains, two unknowns nonlinear variational problem

202
views
0
6 months ago by
Hi FEniCS Community,

I'm coming to you once again with an issue in my problem that I am unable to solve by myself. I am working on the implementation of the two temperature model in a water / gold system. An incoming laser pulse heats up the gold surface, and I want to calculate the values of the electronic and lattice temperatures $T_e$ and $T_l$, respectively, from 0 to 100 ps. A little more details, including the equations, are available in this post.

I have now implemented the suggestions made by Minas Mattheakis found in this other question. However, running my code, I get the following 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

Would you have suggestions about how to progress on this issue? Many thanks for your time, it is appreciated.

The code is:

### Implementing the two temperature model (TTM)

from fenics import *
import numpy as np

### Our special mesh...
def mesh1(lx, ly, Nx, Ny, k):
m = UnitSquareMesh(Nx, Ny)
x = m.coordinates()

x[:,1] = np.sinh(k*(x[:,1] - .5))
x[:,1] -= x[:,1].min()
x[:,1] /= x[:,1].max()

#Scale
x[:,0] = x[:,0]*lx
x[:,1] = x[:,1]*ly

return m

### Main problem

set_log_level(DEBUG)
debug1 = False

# 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

C_water = 4.18       # heat capacity, water / J g-1 K-1
K_water = .0059      # thermal conductivity, water / W cm-1 K-1
rho_water = .997     # density, water / g cm-3
Cl_water = C_water * rho_water # J cm-3 K-1

# 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_metal  = C_metal * rho_metal # J cm-3 K-1

# 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
pulse_diam = .1      # pulse diameter / cm

# Create mesh and define function space
nx = 16
ny = 128
width = 1e-4
height = 1e-4
mesh = mesh1(width, height, nx, ny, 20)

### Define Function Space
P1 = FiniteElement('CG', triangle, 1)
# P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

### Define Functions
u_D = Constant((initial_temperature, initial_temperature))

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

# Define initial value
u_n = interpolate(Constant((302, 300)), V)

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

### Source term:
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=pulse_diam/2., 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 boundary condition
tol = 1E-6 * height
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.sub(0), initial_temperature, boundary_markers, 0))
bcs.append(DirichletBC(V.sub(1), initial_temperature, 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)

if debug1:

# Print all vertices that belong to the boundary parts
for x in mesh.coordinates():
if boundary.inside(x, True):
print('Outer boundary')
print('%s is on boundary' % x)

# Print the Dirichlet conditions
print('Number of Dirichlet conditions:', len(bcs))
if V.ufl_element().degree() == 1:  # P1 elements
d2v = dof_to_vertex_map(V)
coor = mesh.coordinates()
for i, bc in enumerate(bcs):
print('Dirichlet condition %d' % i)
boundary_values = bc.get_boundary_values()
for dof in boundary_values:
print('   dof %2d: u = %g' % (dof, boundary_values[dof]))
if V.ufl_element().degree() == 1:
print('    at point %s' %
(str(tuple(coor[d2v[dof]].tolist()))))

class ThermCoeff(Expression):
def __init__(self, materials, c_0, c_1, **kwargs):
self.materials = materials
self.c_0 = c_0
self.c_1 = c_1

def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.c_0
else:
values[0] = self.c_1

class SourceCoeff(Expression):
"""True if domain is 0"""
def __init__(self, materials, **kwargs):
self.materials = materials

def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = 1
else:
values[0] = 0

Ce_s = ThermCoeff(materials, Ce_, 0, degree=0)
Cl_s = ThermCoeff(materials, Cl_metal, Cl_water, degree=0)
k = SourceCoeff(materials, degree=0)

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

### Write down the equations

Fe0 = Ce_s * Te_n * Te_n * ve * dx(0) \
+ dt * k * f * ve * dx(0) \
- Ce_s * Te * Te * ve * dx(0) \
- dt * k * keq * Te / Tl * inner(grad(Te), grad(ve)) * dx(0) \
- dt * k * G * (Te - Tl) * ve * dx(0)

Fl0 =  Cl_s * Tl_n * vl * dx(0) \
- Cl_s * Tl * vl * dx(0) \
+ dt * k * G * (Te - Tl) * vl * dx(0)

Fe1 = Ce_s * Te_n * Te_n * ve * dx(1) \
+ dt * k * f * ve * dx(1) \
- Ce_s * Te * Te * ve * dx(1) \
- dt * k * keq * Te / Tl * inner(grad(Te), grad(ve)) * dx(1) \
- dt * k * G * (Te - Tl) * ve * dx(1)

Fl1 =  Cl_s * Tl_n * vl * dx(1) \
- Cl_s * Tl * vl * dx(1) \
+ dt * k * G * (Te - Tl) * vl * dx(1)

### Time stepping
steps = np.logspace(-12.2, -10, num_steps+1)
t = 1e-15
dt = steps[1] - steps[0]
i = 0
path = "/home/fenics/shared/TTM results/ttm_test"
res_file = File(path + ".pvd")
res_file_e = File(path + "_e.pvd")
res_file_l = File(path + "_l.pvd")

for t in steps[:-1]:
print('time:', t)
f.t = t
print("  source:", f(Point(width/2., height/2.-1e-9)))

f.t = t

F = action(Fe0 + Fl0 + Fe1 + Fl1, u_n)
J = derivative(F, u_n, u)

## Compute solution
problem = NonlinearVariationalProblem(F, u_n, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

## Compute error
u_e = interpolate(u_D, V)
error = np.abs(u_e.vector().array() - u_n.vector().array()).max()
print('  error: %.3g' % error)

## Update previous solution ???
#     u_n.assign(u)
#     Te_n.assign(Te)
#     Tl_n.assign(Tl)

## Save solution in VTK format
res_file << (u_n, t)
rTe, rTl = u_n.split(True)
res_file_e << (rTe, t)
res_file_l << (rTl, t)

## Increment the time for next step
t += dt
dt = steps[i+1] - steps[i]
i += 1
​

And the debug output (last 2 Newton iterations) gives:
Newton iteration 48: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
 System assembler: boundary condition 0 applies to axis 0
 Elapsed wall, usr, sys time: 1.6492e-05, 0, 0 (DirichletBC init facets)
 Elapsed wall, usr, sys time: 0.000387407, 0, 0 (DirichletBC compute bc)
 System assembler: boundary condition 1 applies to axis 0
 Elapsed wall, usr, sys time: 1.4193e-05, 0, 0 (DirichletBC init facets)
 Elapsed wall, usr, sys time: 0.000344428, 0, 0 (DirichletBC compute bc)
 Elapsed wall, usr, sys time: 9.1255e-05, 0, 0 (Apply (PETScMatrix))
 Elapsed wall, usr, sys time: 0.0940586, 0.1, 0 (Assemble system)
 NewtonSolver: using Jacobian as preconditioner matrix
 Elapsed wall, usr, sys time: 1.5393e-05, 0, 0 (Apply (PETScVector))
 Solving linear system of size 4386 x 4386 (PETSc LU solver, umfpack).
 PETSc Krylov solver starting to solve 4386 x 4386 system.
 Elapsed wall, usr, sys time: 0.0185593, 0, 0.01 (PETSc Krylov solver)
 Elapsed wall, usr, sys time: 0.0187595, 0.01, 0.01 (LU solver)
 Elapsed wall, usr, sys time: 1.3694e-05, 0, 0 (Apply (PETScVector))
 System assembler: boundary condition 0 applies to axis 0
 Elapsed wall, usr, sys time: 1.5692e-05, 0, 0 (DirichletBC init facets)
 Elapsed wall, usr, sys time: 0.00034003, 0, 0 (DirichletBC compute bc)
 System assembler: boundary condition 1 applies to axis 0
 Elapsed wall, usr, sys time: 1.5892e-05, 0, 0 (DirichletBC init facets)
 Elapsed wall, usr, sys time: 0.000309545, 0, 0 (DirichletBC compute bc)
 Elapsed wall, usr, sys time: 3.1685e-05, 0, 0 (Apply (PETScVector))
 Elapsed wall, usr, sys time: 0.106959, 0.1, 0 (Assemble system)
 Newton iteration 49: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
 System assembler: boundary condition 0 applies to axis 0
 Elapsed wall, usr, sys time: 1.949e-05, 0, 0 (DirichletBC init facets)
 Elapsed wall, usr, sys time: 0.000525237, 0, 0 (DirichletBC compute bc)
 System assembler: boundary condition 1 applies to axis 0
 Elapsed wall, usr, sys time: 2.4488e-05, 0, 0 (DirichletBC init facets)
 Elapsed wall, usr, sys time: 0.000538131, 0, 0 (DirichletBC compute bc)
 Elapsed wall, usr, sys time: 0.000106947, 0, 0 (Apply (PETScMatrix))
 Elapsed wall, usr, sys time: 0.0987402, 0.1, 0 (Assemble system)
 NewtonSolver: using Jacobian as preconditioner matrix
 Elapsed wall, usr, sys time: 1.4893e-05, 0, 0 (Apply (PETScVector))
 Solving linear system of size 4386 x 4386 (PETSc LU solver, umfpack).
 PETSc Krylov solver starting to solve 4386 x 4386 system.
 Elapsed wall, usr, sys time: 0.0171559, 0.02, 0 (PETSc Krylov solver)
 Elapsed wall, usr, sys time: 0.0173572, 0.02, 0 (LU solver)
 Elapsed wall, usr, sys time: 1.3293e-05, 0, 0 (Apply (PETScVector))
 System assembler: boundary condition 0 applies to axis 0
 Elapsed wall, usr, sys time: 1.3093e-05, 0, 0 (DirichletBC init facets)
 Elapsed wall, usr, sys time: 0.000313243, 0, 0 (DirichletBC compute bc)
 System assembler: boundary condition 1 applies to axis 0
 Elapsed wall, usr, sys time: 5.697e-06, 0, 0 (DirichletBC init facets)
 Elapsed wall, usr, sys time: 0.000277861, 0, 0 (DirichletBC compute bc)
 Elapsed wall, usr, sys time: 3.2184e-05, 0, 0 (Apply (PETScVector))
 Elapsed wall, usr, sys time: 0.104888, 0.1, 0 (Assemble system)
 Newton iteration 50: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Community: FEniCS Project