What changed between using the linear solver and the nonlinear solver?

General
1

New to Fenics and to numerical methods for PDEs in general. Any insight would be greatly appreciated.

I’m trying to solve a nonlinear problem with mixed boundary conditions (reaction diffusion with Dirichlet, Neumann, and Robin boundary conditions). To start, I thought I would try just solving the Poisson equation using a nonlinear solver, and build everything up from there.

The following code works just fine:

from fenics import *


box_length = 1.0
#finite element parameters
box_length_elements = 128


# Define mesh
mesh = RectangleMesh(Point(0,0), Point(box_length,box_length), box_length_elements, box_length_elements, "right/left")



# Create classes for boundaries
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], box_length)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], box_length)
left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, 1)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

# Define model parameters
D = Constant(1.0)
g_L = Expression("- 10*exp(- pow(x[1] - 0.5, 2))", degree=2)
g_R = Constant(6.0)
f = Constant(1.0)

# Define function space
V = FunctionSpace(mesh, "CG", 2)
u = TrialFunction(V)
v = TestFunction(V)

# Dirichlet boundary conditions
bcs = [DirichletBC(V, 5.0, boundaries, 2),
       DirichletBC(V, 0.0, boundaries, 4)]

# Define measures for boundary terms
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Define variational form
F = inner(D*grad(u), grad(v))*dx - g_L*v*ds(1) - g_R*v*ds(3) - f*v*dx
a, L = lhs(F), rhs(F)

# Solve problem
u = Function(V)
solve(a == L, u, bcs)

# Plot solution and gradient
plot(u)

When I change a == L to F == 0 and correspondingly change u = TrialFunction(V) to u = Function(V), DOLPHIN returns the error

Unable to successfully call PETSc function 'MatSetValuesLocal'

What is going on here? What’s the minimal change to this code to use the nonlinear solver?

2

You probably redefined u in your code.
Here is a minimal example uses the non-linear solver, where the only addition I’ve made to your code is setting a0=1, as you did not define it in your code.

from fenics import *


box_length = 1.0
# finite element parameters
box_length_elements = 128


# Define mesh
mesh = RectangleMesh(Point(0, 0), Point(box_length, box_length),
                     box_length_elements, box_length_elements, "right/left")


# Create classes for boundaries
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], box_length)


class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)


class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], box_length)


left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, 1)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
a0 = 1
# Define model parameters
D = Constant(1.0)
g_L = Expression("- 10*exp(- pow(x[1] - 0.5, 2))", degree=2)
g_R = Constant(6.0)
f = Constant(1.0)

# Define function space
V = FunctionSpace(mesh, "CG", 2)
u = Function(V)
v = TestFunction(V)

# Dirichlet boundary conditions
bcs = [DirichletBC(V, 5.0, boundaries, 2),
       DirichletBC(V, 0.0, boundaries, 4)]

# Define measures for boundary terms
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Define variational form
F = inner(a0*grad(u), grad(v))*dx - g_L*v*ds(1) - g_R*v*ds(3) - f*v*dx
#a, L = lhs(F), rhs(F)

# Solve problem
# u = Function(V)
# solve(a == L, u, bcs)
solve(F == 0, u, bcs=bcs)

# Plot solution and gradient
plot(u)
3

Thanks so much for the quick reply! Yes, it looks like redefining u was the issue. Thanks so much!