NonlinearVariationalSolver converges with zero iteration
I have used NonlinearVariationalSolver as the solver, but what I get is the output "Out: (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
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
- Adjusting the tolerance
- 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, tol) or near(x, width, tol) or near(x, 0, tol) or near(x, 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 <= height / 2. + tol class Omega_1(SubDomain): def inside(self, x, on_boundary): return x >= 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-bw)/w, 2)) * exp(-pow(((t-t0)/tau), 2)) * exp(-alpha*(-x+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()
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.