Taking Derivatives

1

Hello everyone,

I am solving a Poisson equation, and that works well, but I also want to take some partial derivatives from the variational form. I am encountering multiple errors and would really appreciate some help. I am taking two derivatives one with respect to K and one with respect to G.

  1. When I use dK = TrialFunction I get

NotImplementedError: Cannot take length of non-vector expression.

It seems to work if I use dK = Function(V), I do not understand why? Similarly for G.

  1. I can assemble the derivative with respect to K but not with respect to G. When I try to assemble with respect to G I get the following error

ufl.algorithms.check_arities.ArityMismatch: Integrand arguments (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 0, None),) differ from form arguments (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 0, None), Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 1, None)).

Is that because G is on the boundary (ds)?

Thank you very much for the help in advance!

def solve_steady_state_heat(V,mesh,boundary_markers, K, G):
    bmesh = dl.BoundaryMesh(mesh, "exterior", True)
    bdim  = bmesh.topology().dim()
    bot_top_mesh = bot_mesh(mesh,boundary_markers,bmesh,bdim)
    boundary_mesh = bmesh #for ease going from 2D to actual mesh
    submesh_bottom = bot_top_mesh[0]
    submesh_top = bot_top_mesh[1]
    # boundary measure and normal vector
    # ds(1) --> bottom, ds(2) --> top
    ds = dl.Measure("ds", domain=mesh, subdomain_data=boundary_markers)
    # define Trial and Test Functions
    u_trial = dl.TrialFunction(V)
    u_test = dl.TestFunction(V)
    # initialize input functions
    f = dl.Constant(0.0) 
    u0 = dl.Constant(253.0) 
    bmesh_markers = dl.MeshFunction('size_t', bmesh, bdim)
    bmesh_markers.set_all(0)
    # 1 -- bottom
    # 2 -- top
    for i, facet in enumerate(dl.entities(bmesh, bdim)):
        p_meshentity = bmesh.entity_map(bdim)[i]
        p_boundarynumber = boundary_markers.array()[p_meshentity]
        bmesh_markers.array()[i] = p_boundarynumber
    bottom_submesh = dl.SubMesh(bmesh, bmesh_markers, 1)
    top_submesh = dl.SubMesh(bmesh, bmesh_markers, 1)
    bc = dl.DirichletBC(V, u0, boundary_markers, 2)
    
    def a(K,G,u_trial,u_test,f):
        return dl.inner(ufl.exp(K) * dl.grad(u_trial), dl.grad(u_test)) * dl.dx - f * u_test * dl.dx - G*u_test*ds(1)
    res_form = steady_state_heat_varf(Kappa,G,u_trial,u_test,f)
#    dK = dl.TestFunction(V) <-- when. I use I get the error
#   dG = dl.TestFunction(V)
    dK = dl.Function(V) #<--- Question 1) Why does this work? I thought we should use            TestFunction
    dG = dl.Function(V)
    theta_kappa_form = dl.derivative( res_form, K,dK)
    theta_kappa = dl.assemble(theta_kappa_form)
    theta_G_form = dl.derivative( res_form, G,dG)
    theta_G = dl.assemble(theta_G_form)
    
    A_form = ufl.lhs(res_form)
    b_form = ufl.rhs(res_form)
    A, b = dl.assemble_system(A_form, b_form, bcs=bc)
    theta = dl.Function(V)
    dl.solve(A, theta.vector(), b)
    
    return theta
2

Please make sure that the code you supply is a standalone code, that when copy pasted by other can reproduce your error message.

Please also remove all unused variables from the code.

2 Likes
3

Thank you for your feedback Prof. Dokken,

I would directly edit the post to improve it, I cannot find an edit option. I modified one of the basic demos to give me the same errors.

import dolfin as dl
import matplotlib.pyplot as plt
import ufl
from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",element=V.ufl_element())
g = Expression("sin(5*x[0])",element=V.ufl_element())
a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds


# Compute solution
K = Function(V)
solve(a == L, K, bc)

def aa(u,v,f,g):
    return dl.inner(ufl.exp(K)*dl.grad(u), dl.grad(v))*dl.dx - f*v*dl.dx + g*v*ds
res_form = aa(u,v,f,g)

#dK = dl.TestFunction(V) #<--- does not work, it seems we have to use Function, but I guess K is a function as well
dK = dl.Function(V)
aa_K = dl.derivative( res_form, K,dK)
Ak = dl.assemble(aa_K)


#dg = dl.TestFunction(V) #Similar issue
dg = dl.Function(V)
aa_g = dl.derivative( res_form, g,dg)
Ag = dl.assemble(aa_g)
4

In your variational form, you already use a TrialFunction (u) and a TestFunction (v). You cannot have two test functions multiplied by each other in the same form, and therefore the derivative command does not work.

2 Likes
5

I understand so we have multiple Trial and Test functions. Then what is the proper way to take the derivative of the weak form aa with respect to g?

6

In FEniCS there is a concept of scalar, vector and tensors (of order 2).
i.e.
Consider a functional, \int_\Omega c u^2 \mathrm{d}x where u and c is a function.
If we want to compute the derivative (Gateaux derivative, Integrating a functional, followed by differentiation - #2 by kamensky)
we call

from dolfin import *
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
J = c * inner(u, u) * dx
dJdu = derivative(J, u)
v = TestFunction(V)
dJdu_v = derivative(J, u, v)
assert(dJdu == dJdu_v)

assembling dJdu will assemble into a vector:

print(type(assemble(dJdu)))

<class 'dolfin.cpp.la.Vector'>

If you again differentiate dJdu, which is equal to \int_\Omega2 c u\cdot v ~\mathrm{d} x you obtain \int_\Omega 2c \cdot du \cdot v ~\mathrm{d} x which assembles into a matrix:


dJdudv = derivative(dJdu, u)
du = TrialFunction(V)
dJdu_vdu = derivative(dJdu_v, u, du)
assert(dJdudv == dJdu_vdu)
print(type(assemble(dJdudv)))

<class 'dolfin.cpp.la.Matrix'>

In dolfin, we do not have a concept of tensors with higher order than two, and thus you cannot differentiate this form with respect to a TrialFunction or TestFunction, as it would increase the dimensionality of the object.

However, you can differentiate it with respect to a Function, i.e derivative(dJdudv, c, Function(V)) will give you a resulting matrix, and you can create the tensor by looping over all dofs, inserting a one for each dof and zero for all the others if you want to emulate differentiating with respect to a trial or test function.

2 Likes