Linear problem with automatic differentiation

variational formulation
1

Hello,

I would like to define a linear problem using automatic differentiation to define the bilinear form directly from the energy. I wrote the following code

Y = 1.e3
nu = 0.3
L = 1
W = 0.1
mu = Y/(2*(1 + nu))
lmbda = Y*nu/((1 + nu)*(1 - 2*nu))

domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, -W/2]), np.array([L, W/2])], [20, 6], cell_type=mesh.CellType.triangle)

V = fem.VectorFunctionSpace(domain, ("CG", 1))

dw = ufl.TestFunction(V)
w = dolfinx.fem.Function(V)

T = fem.Constant(domain, ScalarType((0, 0)))
ds = ufl.Measure("ds", domain=domain)
dx = ufl.Measure("dx", domain=domain)

def epsilon(u):
    return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(u.geometric_dimension()) + 2*mu*epsilon(u)

energy_density = mu * ufl.inner(epsilon(w), epsilon(w)) + lmbda / 2 * (ufl.tr(epsilon(w))) ** 2

potential_energy = energy_density * dx

bilinear_form = ufl.derivative(potential_energy, w, dw)

linear_form = ufl.inner(T, dw) * dx

fem.petsc.LinearProblem(bilinear_form, linear_form, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

This yields an error:

 117 """Create a PETSc matrix that is compaible with a bilinear form.
    118 
    119 Args:
   (...)
    125 
    126 """
    127 if mat_type is None:
--> 128     return _cpp.fem.petsc.create_matrix(a)
    129 else:
    130     return _cpp.fem.petsc.create_matrix(a, mat_type)

RuntimeError: Cannot create sparsity pattern. Form is not a bilinear form

This is probably because my bilinear is not defined in terms of a ufl.TestFunction and a ufl.TrialFunction, and I don’t know to get such a bilinear form from automatic differntiation.

Many thanks in advance!

2

The bilinear form is actually the second derivative of the energy. You can get this by simply differentiating what is currently (incorrectly) defined as bilinear_form again with respect to the solution w, in the direction of an arbitrary ufl.TrialFunction. (Because this first derivative happens to be linear in w, the second derivative will just replace w with the trial function.)

1 Like
3

thanks! Here is the corrected version of the code

Y = 1.e3
nu = 0.3
L = 1
W = 0.1
mu = Y/(2*(1 + nu))
lmbda = Y*nu/((1 + nu)*(1 - 2*nu))

domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, -W/2]), np.array([L, W/2])], [20, 6], cell_type=mesh.CellType.triangle)

V = fem.VectorFunctionSpace(domain, ("CG", 1))

dw = ufl.TestFunction(V)
hw = ufl.TrialFunction(V)
w = dolfinx.fem.Function(V)

T = fem.Constant(domain, ScalarType((0, 0)))
ds = ufl.Measure("ds", domain=domain)
dx = ufl.Measure("dx", domain=domain)

def epsilon(u):
    return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(u.geometric_dimension()) + 2*mu*epsilon(u)

energy_density = mu * ufl.inner(epsilon(w), epsilon(w)) + lmbda / 2 * (ufl.tr(epsilon(w))) ** 2

potential_energy = energy_density * dx

bilinear_form = ufl.derivative(ufl.derivative(potential_energy, w, dw), w, hw)

linear_form = ufl.inner(T, dw) * dx

problem = fem.petsc.LinearProblem(bilinear_form, linear_form, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})