Manually programmed Newton method doesn't work?

Question

Hi, I have read the following thread Newton method programmed manually and I wanted to try it for my DG method with upwinding for NS equation.

Here is my code

from dolfin import *
import numpy as np

dim = 2

# Mesh
Right   = 5.0
Top = 1.0
mesh    = RectangleMesh(0.0,0.0,Right,Top, 50, 10,"crossed")
bndry   = FacetFunction("size_t", mesh)

# Boundary description
class LeftBoundary(SubDomain):
    def inside(self, x, bndry):
        return x[0] < DOLFIN_EPS

class RightBoundary(SubDomain):
    def inside(self, x, bndry):
        return x[0] > Right-DOLFIN_EPS

class TopBoundary(SubDomain):
    def inside(self, x, bndry):
        return x[1] > Top-DOLFIN_EPS

class BottomBoundary(SubDomain):
    def inside(self, x, bndry):
        return x[1] < DOLFIN_EPS

# Function spaces
V = VectorFunctionSpace(mesh, 'DG', 1)
Q = FunctionSpace(mesh, 'DG', 0)
W = MixedFunctionSpace([V,Q])

dirichlet   = Expression(("0.0", "0.0"), t = 0.0)
u_IN        = Expression(("x[1]*(1.0-x[1])", "0.0"), t = 0.0)

bc_top      = DirichletBC(W.sub(0), dirichlet, TopBoundary(), 'pointwise')
bc_bottom   = DirichletBC(W.sub(0), dirichlet, BottomBoundary(), 'pointwise')
bc_left     = DirichletBC(W.sub(0), u_IN, LeftBoundary(), 'pointwise')

bcs = [bc_top, bc_bottom, bc_left]

# Unknown and test functions
w   = Function(W)
(u, p)  = split(w)
(v, q)  = TestFunctions(W)
u_1 = Function(V) # velocity from previous time step


# Problem parameters
T_end       = 1.0
dt      = 0.1
u0      = Expression(("0.0", "0.0"))
f_rhs       = Expression(("0.0", "0.0"), t = 0.0)

sigma   = 10.0

# Variational formulation
n   = FacetNormal(mesh)
ds  = Measure("ds", subdomain_data = bndry)
I   = Identity(u.geometric_dimension())
D   = 0.5*(grad(u) + grad(u).T)
def F(D):
    F = D
    return F

def a(u,v):
    M = inner(F(D),grad(v))*dx - inner(avg(F(D))*n('+'),jump(v))*dS 
    return M

def J(u,v):
    M = sigma*inner(jump(u),jump(v))*dS
    return M

def b(p,v):
    M = -p*div(v)*dx  + avg(p)*dot(jump(v),n('+'))*dS
    return M

def c(u,w,v):   
    P = avg(dot(w,n))
    H = conditional(P < 0.0, dot(u('+'),w('+')), dot(u('-'),w('-')))
    M = -0.5*inner(grad(v)*u,w)*dx + inner(0.5*H*n('+'),jump(v))*dS
    return M

def L(v):
    M = inner(f_rhs,v)*dx 
    return M

T =  (1/dt)*inner(u - u_1,v)*dx + a(u,v) + J(u,v) + b(p,v) + b(q,u) + c(u,u_1,v) - L(v) 

dup = TrialFunction(W)
Jac     = derivative(T, w, dup)

# Solution

plt = plot(u, mesh = mesh, mode="color", interactive = True, wireframe = False)
pltP = plot(p, mesh = mesh, mode="color", interactive = True, wireframe = False)


a_tol, r_tol = 1e-7, 1e-10
w_inc = Function(W)
u_1.assign(interpolate(u0,V))

bc_top_du   = DirichletBC(W.sub(0), Constant((0.0,0.0)), TopBoundary(), 'pointwise')
bc_bottom_du    = DirichletBC(W.sub(0), Constant((0.0,0.0)), BottomBoundary(), 'pointwise')
bc_left_du  = DirichletBC(W.sub(0), Constant((0.0,0.0)), LeftBoundary(), 'pointwise')

bcs_du = [bc_top_du, bc_bottom_du, bc_left_du]

ufile = File("u.xdmf")
pfile = File("p.xdmf")
t = dt
u_1.rename("u", "velocity")
ufile << u_1

while t<=T_end:
    dirichlet.t     = t
    u_IN.t      = t
    nIter = 0
    eps = 1
    print 'time step: ', t
    if True:
        problem = NonlinearVariationalProblem(T, w, bcs, Jac)
        solver = NonlinearVariationalSolver(problem)
        solver.solve()
    else:
        while eps > 1e-10 and nIter < 4:    # Newton iterations
        nIter += 1
            A, b = assemble_system(Jac, -T, bcs_du)
            solve(A, w_inc.vector(),b)     # Determine step direction
            eps = np.linalg.norm(w_inc.vector().array(), ord = 2)
            a = assemble(T)
            for bc in bcs_du:
                bc.apply(a)
            print 'b.norm("l2") = ', b.norm('l2'), 'np.linalg.norm(a.array(), ord = 2) = ', np.linalg.norm(a.array(), ord = 2)
            fnorm = b.norm('l2')
            lmbda = 1.0     # step size, initially 1
            w.vector()[:] += lmbda*w_inc.vector()    # New w vector
        (u,p) = w.split(True)
        print '      {0:2d}  {1:3.2E}  {2:5e}'.format(nIter, eps, fnorm)    

    (u,p) = w.split(True)
    u_1.assign(u)       
    u.rename("u", "velocity")
    p.rename("p", "pressure")
    ufile << (u,t)
    pfile << (p,t)
    t += float(dt)
plt.plot(u)
pltP.plot(p)

If you leave "True" in the if statement, regular FENICS' Newton is called and you should get proper solution (nonzero solution for given data and problem). If you rewrite it to "False", manually writen Newton is performed.
The problem is that this time You get zero solution.

Can you pls help me find the right way to implement Newton here? I.e. what should I change in my implementation to get the same results as I get with FENICS' Newton solver for the same problem?

Comments

It is not clear what you want to ask. From what I understand your Newton implementation gives the right result for homogeneous Dirichlet conditions and fails to converge to the right solution for inhomogeneous ones?

---

Sry, I hope it is now more clear. I just wonder how should I implement Newton to work properly for my problem.

Answer 1

The mistake seems to be that in your "hand-made" implementation, you don't apply your boundary condition ("bcs" in your code) to the initial solution w. So you start your Newton iteration with w = 0 and then your Newton increments also have a zero boundary condition (bcs_du), so there is no way that a nonzero boundary condition can enter your equation. This is why you get the correct solution, which is zero in this case.
What you have to do is something like

for bc in bcs:
    bc.apply(w)

after you define w at the beginning.
By the way, posting 150 lines of code is not really considered good practice ;)
Cheers, Gregor

Comments

Thank you very much. Everything works fine and is clear to me. It just suffices to add

for bc in bcs:
    bc.apply(w.vector())

as You proposed. Next time, I will try to make the code as short as possible. ;)