Hello. I have been trying to solve a more complex problem, but decided to first prove my code works for a 1D time-independent problem that is non-linearly coupled between two unknown functions, phi and c. The two equations are
-\kappa \partial ^2 \phi / \partial x^2 = \gamma * c* exp(\alpha (\phi_1 - \phi - \phi_0 )
D d ^2 c / dx^2 = - \beta * d ^2 \phi /d x^2
with boundary conditions on domain x [0, L]
c(0) = c_0, \phi (0) = \phi_0, dc(L)/dx = 0, d phi[L] /dx = 0
kappa, gamma, D, beta, phi1, and phi0 are constants.
Combining these equations and multiplying the first by test function v_1 and second by test function v_2, the variational problem becomes
\int (\beta dot(\nabla \phi , \nabla v_1) + D dot( \nabla c, \nabla v_1))dx +
\int (\kappa dot(\nabla \phi , \nabla v_2)dx - \gamma * c * exp(\alpha * (\phi_1 - \phi - \phi_0) ) v_2 )dx
However, I get this 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
***
*** DOLFIN version: 2019.1.0
*** Git changeset: unknown
I’m running this on Ubuntu 20.04 LTS (GNU/Linux 4.4.0-19041-Microsoft x86_64). My code is attached:
from dolfin import *
import mshr
L = .1 # Length of electrode, cm
kap = .292 # Corrected electrolyte conductivity, S/cm
n = 1 # Electrons transfered
alpha = 19.5 # exponent coefficient
D = 5 *10**-5 # Diffusion coefficient, cm^2/s
Far = 9485 # Faraday constant, C/mol e
ai0 = 2.5 # exchange current density, A/cm^3
c_0 = 1 # initial bulk concentration, mol
phi0 = 0 # equilibrium potential, V
phi1 = .5 # applied potential, V
# Constants in current relation
gamma = Constant(ai0/c_0) # coefficient
beta = Constant(kap/(n*Far)) # current coeff in continuity
# Define boundaries; made 1D, getting rid of y-direction
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], L)
# Initialize sub-domain instances
left = Left()
right = Right()
# Create mesh.
mesh = IntervalMesh(100, 0, L)
# Marking boundary dCmains.
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
CG1_elem = FiniteElement("CG", mesh.ufl_cell(), 1)
#Defining the function spaces for the concentration
V_c = FunctionSpace(mesh, CG1_elem)
#Defining the function spaces for the voltage
V_phi = FunctionSpace(mesh, CG1_elem)
#Defining the mixed function space
W_elem = MixedElement([CG1_elem, CG1_elem])
W = FunctionSpace(mesh, W_elem)
#Defining the "Trial" functions
u = Function(W)
c, phi = split(u)
#Defining the test functions
(v_1, v_2) = TestFunctions(W)
# Newman condition automatically applied at right boundary.
# Define Dirichlet boundary conditions at left boundary.
bc_phi = DirichletBC(W.sub(1), phi0, boundaries, 1)
bc_con = DirichletBC(W.sub(0), c_0, boundaries, 1)
#Collecting the boundary conditions
bcs = [bc_phi, bc_con]
# Define variational form
F = c * dot(grad(phi), grad(v_2)) * dx \
- gamma*exp(alpha*(phi1 - phi - phi0)) * c * v_2 * dx \
+ D * dot(grad(c),grad(v_1)) * dx \
+ beta * dot(grad(phi),grad(v_1)) * dx
# Create VTK files for visualization output
vtkfile_phi = File('ss_1D/phi.pvd')
vtkfile_c = File('ss_1D/c.pvd')
solve(F == 0, u, bcs)
# Save solution to file (VTK)
vtkfile_c << (c,)
vtkfile_phi << (phi)
I’ve noticed other people post with similar convergence issues and am wondering if anyone has made progress. I’ve thought about changing the tolerance of my solver, but it would have to be pretty high, as I’m getting this result when solving:
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 7.217e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 7.014e+02 (tol = 1.000e-10) r (rel) = 9.719e+01 (tol = 1.000e-09)
Newton iteration 2: r (abs) = inf (tol = 1.000e-10) r (rel) = inf (tol = 1.000e-09)
Newton iteration 3: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)