Help with Gateaux derivative

General
1

I am working on a non-linear heat equation with the following bilinear form:

a(u^{n+1}, v) = \int c(u^{n+1}) \frac{u^{n+1} - u^n}{\Delta t} v + k \int \nabla u^{n+1} \cdot \nabla v = 0

where u^{n+1} is the solution we seek and v is the test function. The residual is just the bilinear form evaluated for a particular u^{n+1} = u:

F(v) = \int c(u) \frac{u - u^n}{\Delta t} v + k \int \nabla u \cdot \nabla v = 0 
F = c(unp1)*(unp1 - un)/dt*v*ufl.dx + k * ufl.dot(ufl.grad(unp1),ufl.grad(v)) * ufl.dx

and I believe the Jacobian is:

J(\delta u, v) = \int c(u^{n+1}) \frac{\delta u}{\Delta t} v + \int c'(u^{n+1}) \frac{\delta u}{\Delta t} u^{n+1} v + k \int \nabla \delta T \cdot \nabla v 
  J  =  c(unp1)/dt* delta_u * v * ufl.dx
  J +=  cp(unp1)/dt * delta_u * unp1 * v * ufl.dx
  J += k * ufl.dot( ufl.grad(delta_u), ufl.grad(v) ) * ufl.dx

However when I plug this Jacobian and I run Newton-Raphson, the result is incorrect / different than if I used ufl.derivative. I know because ufl.derivative converges much faster and other examples with non-linearity in the conductivity showed no difference between the Jacobian I wrote and the one of ufl.derivative. I suspect this line is wrong but I can’t see why:

J +=  cp(unp1)/dt * delta_u * unp1 * v * ufl.dx

Could you help me write the right Jacobian form ?

MWE:

import dolfinx
import dolfinx.fem.petsc
import matplotlib.pyplot as plt
import numpy as np
import ufl
from mpi4py import MPI
from petsc4py import PETSc
import matplotlib.pyplot as plt

N = 10
num_timesteps = 1

domain = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, N)
V = dolfinx.fem.functionspace(domain, ("Lagrange", 1))
unp1  = dolfinx.fem.Function(V)
un    = dolfinx.fem.Function(V)

# IC
f0 = lambda x : 2 - 0.5*(np.cos(3*x[0]) - np.sin(5*x[0]))
unp1.interpolate(f0)

dt = dolfinx.fem.Constant(domain,0.1)
k = 1.0
c = lambda T : 2.0 + T**2
cp = lambda T : 2*T

v = ufl.TestFunction(V)
F = c(unp1)*(unp1 - un)/dt*v*ufl.dx + k * ufl.dot(ufl.grad(unp1),ufl.grad(v)) * ufl.dx
residual = dolfinx.fem.form(F)

own_derivative = True
J = ufl.derivative(F, unp1)
if own_derivative:
    delta_u = ufl.TrialFunction(V)
    J  =  c(unp1)/dt* delta_u * v * ufl.dx
    J +=  cp(unp1)/dt * delta_u * unp1 * v * ufl.dx
    J += k * ufl.dot( ufl.grad(delta_u), ufl.grad(v) ) * ufl.dx
jacobian = dolfinx.fem.form(J)

A = dolfinx.fem.petsc.create_matrix(jacobian)
L = dolfinx.fem.petsc.create_vector(residual)

solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
du = dolfinx.fem.Function(V)

coords = V.tabulate_dof_coordinates()[:, 0]
sort_order = np.argsort(coords)
max_iterations = 25
solutions = np.zeros((max_iterations + 1, len(coords)))
solutions[0] = unp1.x.array[sort_order]

for time_iter in range(num_timesteps):
    i = 0
    un.x.array[:] = unp1.x.array[:]
    while i < max_iterations:
        # Assemble Jacobian and residual
        with L.localForm() as loc_L:
            loc_L.set(0)
        A.zeroEntries()
        dolfinx.fem.petsc.assemble_matrix(A, jacobian)
        A.assemble()
        dolfinx.fem.petsc.assemble_vector(L, residual)
        L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)

        # Scale residual by -1
        L.scale(-1)
        L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

        # Solve linear problem
        solver.solve(L, du.x.petsc_vec)
        du.x.scatter_forward()
        # Update u_{i+1} = u_i + delta u_i
        unp1.x.array[:] += du.x.array
        i += 1

        # Compute norm of update
        correction_norm = du.x.petsc_vec.norm(0)
        print(f"Iteration {i}: Correction norm {correction_norm}")

        solutions[i, :] = unp1.x.array[sort_order]
        if correction_norm < 1e-10:
            break
2

I think you missed a term in your derivative:

J(\delta u, v) = \int c(u^{n+1}) \frac{\delta u}{\Delta t} v \; \mathrm{d}x + \int c'(u^{n+1}) \frac{\delta u}{\Delta t} \left(u^{n+1} - u^n\right) v \; \mathrm{d} x + k \int \nabla \delta u \cdot \nabla v \; \mathrm{d} x

I haven’t tried your code, but hopefully this gives quadratic convergence.

2 Likes