Printing out the Jaccobian and residual

1

Dear colleagues,
In my code I am dealing with non linear problem and would like to print out the Jaccobain matrix and residual vector from variational form. When I try using the assemble_system it return the following error:

*** Error:   Unable to assemble system.
*** Reason:  expected a linear form for L.
*** Where:   This error was encountered inside SystemAssembler.cpp.
*** Process: 0

Any alternative suggestion or correction would be appreciated if my approach is altogether wrong.
My simplified problem code is following:


u_D = Constant(100)

mesh = UnitCubeMesh(5,5,5)

V = FunctionSpace(mesh, 'P', 1)

# Define initial value
u_n = interpolate(u_D, V)

v = TestFunction(V)
f = Constant(0)


# List
boundary_conditions = {2: {'Dirichlet':  u_D}}

class BoundaryZ0(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0, 1e-12)


# Mark boundaries
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 9999)
by0 = BoundaryZ0(); by0.mark(boundary_markers, 2)

# Collect Dirichlet conditions
bcs = []
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc_cur = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
                             boundary_markers, i)
        bcs.append(bc_cur)


u = Function(V)
dt = 1e-2

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx 


u.interpolate(u_D)

# Time-stepping
t = 0

for n in range(2):
    t += dt

    J = derivative(F, u);

    solve(F == 0, u, bcs)

    if n == 0:
        # Solver parameters
        prm=parameters
        info(prm, True)

    # Update previous solution
    u_n.assign(u)
    F_u_n = action(F,u_n);  # F(u_0)
    A, b = assemble_system(J, F_u_n, bcs)
2

I do not see the nonlinearity. Here is what I believe to be a quick fix, that is you need to have test and trial functions to assemble the desired system matrix.

rom dolfin import *

u_D = Constant(100)

mesh = UnitCubeMesh(5,5,5)

V = FunctionSpace(mesh, 'P', 1)

# Define initial value
u_n = interpolate(u_D, V)

v = TestFunction(V)
f = Constant(0)


# List
boundary_conditions = {2: {'Dirichlet':  u_D}}

class BoundaryZ0(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0, 1e-12)


# Mark boundaries
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 9999)
by0 = BoundaryZ0(); by0.mark(boundary_markers, 2)

# Collect Dirichlet conditions
bcs = []
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc_cur = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
                             boundary_markers, i)
        bcs.append(bc_cur)


u = Function(V)
utrial = TrialFunction(V)
dt = 1e-2


a = utrial*v*dx + dt*dot(grad(utrial), grad(v))*dx
L = (u_n + dt*f)*v*dx

u.interpolate(u_D)

# Time-stepping
t = 0

for n in range(2):
    t += dt

    solve(a == L, u, bcs)

    if n == 0:
        # Solver parameters
        prm=parameters
        info(prm, True)

    # Update previous solution
    u_n.assign(u)
    A, b = assemble_system(a, L, bcs)
3

Thanks @tucker ! It worked for me as a quick solution. However, in actual code I have a term depending upon u is involved in variational form which I removed to keep minimum working example simple. Any suggestions if I keep it considering the non-linear problem in form of F=0 and want to print the Jaccobian and residual (because I think I wont be able to access A, b matrices from non-linear problem then) ? Thanks :slight_smile:

4

Try this. Be careful as the Newton system is A x = -1 * b

from dolfin import *

u_D = Constant(100)

mesh = UnitCubeMesh(5,5,5)

V = FunctionSpace(mesh, 'P', 1)

# Define initial value
u_n = interpolate(u_D, V)

v = TestFunction(V)
f = Constant(0)


# List
boundary_conditions = {2: {'Dirichlet':  u_D}}

class BoundaryZ0(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0, 1e-12)


# Mark boundaries
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 9999)
by0 = BoundaryZ0(); by0.mark(boundary_markers, 2)

# Collect Dirichlet conditions
bcs = []
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc_cur = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
                             boundary_markers, i)
        bcs.append(bc_cur)


u = Function(V)
utrial = TrialFunction(V)
dt = 1e-2


F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
u.interpolate(u_D)

# Time-stepping
t = 0

for n in range(2):
    t += dt
    J = derivative(F, u, utrial)
    solve(F == 0, u, bcs)

    if n == 0:
        # Solver parameters
        prm=parameters
        info(prm, True)

    # Update previous solution
    u_n.assign(u)
    A, b = assemble_system(J, F, bcs)