Different Methods for Solving Linear System

General
1

Hi everyone,

I was just going through a tutorial for solving the nonlinear Poisson equation and I noticed something odd. After defining the linear and nonlinear forms as well as the boundary conditions, calling ‘solve(a == L, u, bc)’ produces the correct output.

However, if I instead explicitly create the matrix and right-hand side vector for the linear algebra problem, the program runs but does not give the correct output.

What’s causing the discrepancy?

Working program:

from fenics import *
import numpy as np

nx = ny = 50
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)


tol = 1E-14
def left_boundary(x, on_boundary):
    return on_boundary and abs(x[0]) < tol

def right_boundary(x, on_boundary):
    return on_boundary and abs(x[0]-1) < tol

Gamma_0 = DirichletBC(V, Constant(0.0), left_boundary)
Gamma_1 = DirichletBC(V, Constant(1.0), right_boundary)
bcs = [Gamma_0, Gamma_1]

def q(u):
    return (1+u)**2


u = TrialFunction(V)
v = TestFunction(V)
u_k = interpolate(Constant(0.0), V)  # previous (known) u
a = inner(q(u_k)*grad(u), grad(v))*dx
f = Constant(0.0)
L = f*v*dx

# Picard iterations
u = Function(V)     # new unknown function
eps = 1.0           # error measure ||u-u_k||
tol = 1.0E-5        # tolerance
iter = 0            # iteration counter
maxiter = 25        # max no of iterations allowed
while eps > tol and iter < maxiter:
    iter += 1
    solve(a == L, u, bcs)
    diff = u.vector().get_local() - u_k.vector().get_local()
    eps = np.linalg.norm(diff, ord=np.Inf)
    print('iter=%d: norm=%g' % (iter, eps))
    u_k.assign(u)   # update for next iteration

Malfunctioning program:

from fenics import *
import numpy as np

nx = ny = 50
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)


tol = 1E-14
def left_boundary(x, on_boundary):
    return on_boundary and abs(x[0]) < tol

def right_boundary(x, on_boundary):
    return on_boundary and abs(x[0]-1) < tol

Gamma_0 = DirichletBC(V, Constant(0.0), left_boundary)
Gamma_1 = DirichletBC(V, Constant(1.0), right_boundary)
bcs = [Gamma_0, Gamma_1]

def q(u):
    return (1+u)**2


u = TrialFunction(V)
v = TestFunction(V)
u_k = interpolate(Constant(0.0), V)  # previous (known) u
a = inner(q(u_k)*grad(u), grad(v))*dx
f = Constant(0.0)
L = f*v*dx

assembler = SystemAssembler(a, L, bcs)
solver = LUSolver("mumps")
A = Matrix()
solver.set_operator(A)
assembler.assemble(A)
b = assemble(L)

# Picard iterations
u = Function(V)     # new unknown function
eps = 1.0           # error measure ||u-u_k||
tol = 1.0E-5        # tolerance
iter = 0            # iteration counter
maxiter = 25        # max no of iterations allowed
while eps > tol and iter < maxiter:
    iter += 1
    solve(A, u.vector(), b)
    diff = u.vector().get_local() - u_k.vector().get_local()
    eps = np.linalg.norm(diff, ord=np.Inf)
    print('iter=%d: norm=%g' % (iter, eps))
    u_k.assign(u)   # update for next iteration
2

I can’t run the code right now, but this looks suspicious. You assemble A using the assembler, which will also apply boundary conditions. Instead, you use the assemble free function for b, which will NOT apply boundary conditions.

2 Likes
3

I don’t think this is the issue. If I replace all of

assembler = SystemAssembler(a, L, bcs)
solver = LUSolver("mumps")
A = Matrix()
solver.set_operator(A)
assembler.assemble(A)
b = assemble(L)

with

[A,b] = assemble_system(a,L,bcs)

I still get incorrect results.

4

That must be put inside the for loop, otherwise you aren’t using the updated u_k.

5

Oh goodness. That was silly of me. Thanks

1 Like