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

202

views

0

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:

And the debug output (last 2 Newton iterations) gives:

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

Please login to add an answer/comment or follow this question.