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
Please login to add an answer/comment or follow this question.

Similar posts:
Search »