Solving a non linear poisson problem, where q(u) is function of grad(u)

General
1

I am trying to solve a non-linear poisson’s equation using Newton-Rapson Method

Explicitly -del((K0 + K1(norm(grad(u))**2)).del(u)) = f, where let f(u) = K0 + K1(norm(grad(u))**2
The textbook has an example in which f(u) is a function of u and is (1+u)**m.
I am trying to manually do the newton raphson loop instead of using the non linear variational in fenics.

Accordinly, I have linearized the equation and i am getting,

Screenshot 2023-06-08 113708

Accordinly I have done the code as,

from fenics import *
import numpy as np
import matplotlib.pyplot as plt  

K0 = 200
K1 = K0/10

Nx = 10
Ny = 10

# u = np.zeros(Nx, Ny)

mesh = UnitSquareMesh(Nx, Ny)    # mesh
V = FunctionSpace(mesh, 'Lagrange', 1) # function space

# define top boundary
def top(x, on_boundary):
    return on_boundary and near(x[1], 1)

def bottom(x, on_boundary):
    return on_boundary and near(x[1], 0)

def left(x, on_boundary):  
    return on_boundary and near(x[0], 0)

def right(x, on_boundary):
    return on_boundary and near(x[0], 1)


# Define boundary condition
g = Constant(0.0)

# Top boundary (index 3) is fixed
bc1 = DirichletBC(V, g, top)
bc2 = DirichletBC(V, g, bottom)
bc3 = DirichletBC(V, g, left)
bc4 = DirichletBC(V, g, right)
bcs = [bc1, bc2, bc3, bc4]

# Define Variational problem

f = Constant(-6.0)

u = TrialFunction(V)
v = TestFunction(V)
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx
# Compute solution
A, b = assemble_system(a, L, bcs)

u_k = Function(V)
U_k = u_k.vector()
solve(A, U_k, b)

# Till here, the code involves getting an initial u from a linear poisson equation. Now, to
# do the Newton raphson,

Gamma_0_du = DirichletBC(V, Constant(0), left)
Gamma_1_du = DirichletBC(V, Constant(0), right)
Gamma_2_du = DirichletBC(V, Constant(0), top)
Gamma_3_du = DirichletBC(V, Constant(0), bottom)

bcs_du = [Gamma_0_du, Gamma_1_du, Gamma_2_du, Gamma_3_du]

def q(u):
    return (K0 + K1*dot(grad(u),grad(u))**2)
# def Dq(u):
    # return m*(1+u)**(m-1)


du = TrialFunction(V) # u = u_k + omega*du
#right = inner((K0 + 2*K1* dot(grad(u_k),grad(du)))*nabla_grad(u_k), nabla_grad(v))*dx
a = inner( q(u_k)*nabla_grad(du), nabla_grad(v))*dx + inner((K0 + 2*K1* dot(grad(u_k),grad(du)))*nabla_grad(u_k), nabla_grad(v))*dx
L = - inner(q(u_k)*nabla_grad(u_k), nabla_grad(v))*dx


du = Function(V)
u = Function(V) # u = u_k + omega*du
omega = 1.0 # relaxation parameter
eps = 1.0
tol = 1.0E-5
iter = 0
maxiter = 25
while eps > tol and iter < maxiter:
    iter += 1
    A, b = assemble_system(a, L, bcs_du)
    solve(A, du.vector(), b)
    eps = np.linalg.norm(du.vector().array(), ord=np.Inf)
    print("Norm:", eps)
    u.vector()[:] = u_k.vector() + omega*du.vector()
    u_k.assign(u)

I am getting error as
ArityMismatch Traceback (most recent call last)
Cell In[4], line 29
27 while eps > tol and iter < maxiter:
28 iter += 1
—> 29 A, b = assemble_system(a, L, bcs_du)
30 solve(A, du.vector(), b)
31 eps = np.linalg.norm(du.vector().array(), ord=np.Inf)
ArityMismatch: Adding expressions with non-matching form arguments () vs (‘v_1’,).

Can you please help me with this problem, and see where I am going wrong

Also in addition, with respect to the example given in textbook,

du = Function(V)
u = Function(V) # u = u_k + omega*du
omega = 1.0 # relaxation parameter
eps = 1.0
tol = 1.0E-5
iter = 0
maxiter = 25
while eps > tol and iter < maxiter:
    iter += 1
    A, b = assemble_system(a, L, bcs_du)
    solve(A, du.vector(), b)
    eps = numpy.linalg.norm(du.vector().array(), ord=numpy.Inf)
    print ’Norm:’, eps
    u.vector()[:] = u_k.vector() + omega*du.vector()
    u_k.assign(u)

It pops an error ‘dolfin.cpp.la.PETScVector’ object has no attribute ‘array’. How to deal with this?