# Derivative function has no attribute 'subs'

**URL:** <https://fenicsproject.discourse.group/t/derivative-function-has-no-attribute-subs/4956>\
**Category:** variational formulation\
**Created:** [January 23, 2021, 1:48am UTC](https://fenicsproject.discourse.group/t/derivative-function-has-no-attribute-subs/4956 "2021-01-23T01:48:04Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![mary\_h](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mary_h/32/1943_2.png) [@mary\_h](https://fenicsproject.discourse.group/u/mary_h)\
**Post date:** [January 23, 2021, 1:48am UTC](https://fenicsproject.discourse.group/t/derivative-function-has-no-attribute-subs/4956/1 "2021-01-23T01:48:05Z")

</div>

Hello all,

I am new to Fenics and I want to implement a gradient descent algorithm to minimize the energy in the hyperelasticity demo. I couldn’t find any readily available solver for gradient descent in Fenics, and I tried to implement the algorithm by myself. However, when I used “subs” to substitute y\_initial in the derivative of energy (F in this code), I got this error:  
`

> AttributeError: ‘Form’ object has no attribute ‘subs’

`  
Can anyone help me with solving this problem?  
Following is my minimal code:

```
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

# Create mesh and define function space
mesh = RectangleMesh(Point(-1,-1),Point(1,1), 50, 50)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

################################################################################
# Define boundary condictions
################################################################################

F11=1
F12=0.0 + 0.1
F21=0.0 - 0.1
F22=1.0  

y_BC = Expression(("F11*x[0]+F12*x[1]", "F21*x[0]+F22*x[1]"), F11=F11,F12=F12,F21=F21,F22=F22,degree=2)

tol = 1E-3
def x0_boundary(x, on_boundary):
    return on_boundary and near(x[0], -1, tol) or near(x[0], 1, tol) or \
                            near(x[1], -1, tol) or near(x[1], 1, tol) 

bcs = [DirichletBC(V, y_BC, x0_boundary)]

################################################################################
# Define functions
################################################################################

dy = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
y = Function(V)
B = Constant((0, 0)) # Body force per unit area
T = Constant((0, 0)) # Traction force on the boundary

y_reference = Expression(("x[0]","x[1]"),degree=3)

################################################################################
# Elasticity model: Stored strain energy density (compressible neo-Hookean model)
################################################################################

E, nu = 10.0, 0.3
mu, blk = Constant(E/(2*(1 + nu))), Constant(E/(3*(1 - 2*nu)))

def W(y): #Stored strain energy density 
    F = grad(y) # Deformation gradient
    C = F.T*F # Right Cauchy-Green tensor
    Ic = tr(C)
    J = det(F)
    return (mu/2)*(Ic - 3) - mu*ln(J) + (blk/2 - mu/3)*(J-1)**2

################################################################################
# Total potential energy and solve
################################################################################

Pi = W(y)*dx -dot(B, (y-y_reference))*dx + dot(T, (y-y_reference))*ds
F = derivative(Pi, y, v) # first variation (derivative at u in the direction of v)
dF = derivative(F, y, dy) # Jacobian of F

####gradient descent solver:
alpha_ = 0.01
i=0   
 
y_initial = project(y_BC, V)
y_new = Function(V)    
y_new = project(y_BC, V)
    
while (np.abs(F.subs(y, y_new) ) >1e-6):
    if i!=0:
        y_initial = y_new
            
    y_new = y_initial - alpha_ * J.subs(y, y_initial)    
    i = i+1

y = y_new

```

Thank you so much for your help!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 23, 2021, 11:42am UTC](https://fenicsproject.discourse.group/t/derivative-function-has-no-attribute-subs/4956/2 "2021-01-23T11:42:30Z")

</div>

Note that F is a variational formulation (symbolical) written in the UFL form language.  
To replace a coefficient in a ufl-form, you need to use `ufl.replace`. i.e.

```auto
from ufl import replace

F_ = replace(F, {y:y_new})

```

Note that this returns a ufl-form, a symbolical representation of `F`. To get a numerical value, you need to assemble it:

```auto
F_assembled = assemble(F_)

```

This returns a `dolfin.Vector` which you can take the norm of:

```auto
F_norm = F_assembled.norm("l2")

```

For your specific problem, this would end up as:

```auto
from ufl import replace
while assemble(replace(F, {y:y_new})).norm("l2") >1e-6:
     # Add code

```

> [@mary\_h](#):
>
> ` y_new = y_initial - alpha_ * J.subs(y, y_initial) `

The same hold for this assignment, you need to use replace and assembly.

Please also note that you are constantly reassigning values in your loop.  
You should use the `assign` function :

```auto
y_new = Function(V)
y_new.assign(project(y_BC, V))
while ...:
      if i!=0:
          y_initial.assign(y_new)
     y_new.vector()[:] = y_initial.vector()[:] - assemble(.....)[:]

```

Note that `J` is not defined outside of `W(y)`.

---

<div class="post-metadata">

**Author:** ![bhaveshshrimali](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bhaveshshrimali/32/124_2.png) [@bhaveshshrimali](https://fenicsproject.discourse.group/u/bhaveshshrimali)\
**Post date:** [January 23, 2021, 5:48pm UTC](https://fenicsproject.discourse.group/t/derivative-function-has-no-attribute-subs/4956/3 "2021-01-23T17:48:13Z")

</div>

Maybe I am skipping some details, or the formulation is meant to be a simplified version of another more involved problem, but the energy should be `mu/2. * (Ic - 2)` for `2D` right? You may also want to take a look at [Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation](https://comet-fenics.readthedocs.io/en/latest/demo/2D_plasticity/vonMises_plasticity.py.html) or a post from the old forum [Newton method programmed manually - FEniCS Q&A](https://fenicsproject.org/qa/536/newton-method-programmed-manually/)

Both of them write the Newton’s method manually, but that should be easily modifiable for gradient descent.

In any case, a simplified version of how I would go about implementing gradient descent would be something like (I have used displacement `u` instead of the deformation mapping `y` as in your case)

```auto
from dolfin import *
from copy import deepcopy

mesh = UnitCubeMesh(3,3,3)
V = VectorFunctionSpace(mesh, "CG", 1)
u = Function(V)
v = TestFunction(V)
du = TrialFunction(V)
B = Constant((0., -0.5, 0.))

left = CompiledSubDomain("near(x[0], 0.)")
bcs = [
    DirichletBC(V, Constant((0., 0., 0.)), left)
]

E, nu = Constant(10.), Constant(0.3)
mu = E/2./(1+nu)
lmbda = 2 * mu * nu / (1. - 2* nu)

F = Identity(len(u)) + grad(u)

Pi = mu/2. * (inner(F, F) - 3.) - mu * ln(det(F)) + lmbda/2. * (det(F) - 1.)**2
Pi *= dx
Pi -= inner(B, u) * dx

residual = derivative(Pi, u, v)
jacobian = derivative(residual, u, du)

A, b = assemble_system(jacobian, -residual, bcs)

res_initial = b.norm("l2")
res = deepcopy(res_initial)

alph = 0.05
max_iters = 5000
iters = 1

while res/res_initial > 1.e-6 and iters < max_iters:
    u.vector()[:] += alph * b[:]
    _, b = assemble_system(jacobian, -residual, bcs)
    res = b.norm("l2")
    print(f"Norm: {res}, iter: {iters}")
    iters += 1
    

```

It would be terribly slow, but then again the problem itself could warrant such a method in some cases.

---

<div class="post-metadata">

**Author:** ![mary\_h](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mary_h/32/1943_2.png) [@mary\_h](https://fenicsproject.discourse.group/u/mary_h)\
**Post date:** [January 23, 2021, 6:27pm UTC](https://fenicsproject.discourse.group/t/derivative-function-has-no-attribute-subs/4956/4 "2021-01-23T18:27:12Z")

</div>

Thank you so much for your very helpful responses!

---

<div class="post-metadata">

**Author:** ![mary\_h](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mary_h/32/1943_2.png) [@mary\_h](https://fenicsproject.discourse.group/u/mary_h)\
**Post date:** [January 23, 2021, 6:39pm UTC](https://fenicsproject.discourse.group/t/derivative-function-has-no-attribute-subs/4956/5 "2021-01-23T18:39:00Z")

</div>

Thanks for catching that. Yes, it should be `mu/2. * (Ic - 2)` for my 2D problem.

---

<div class="post-metadata">

**Author:** ![mary\_h](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mary_h/32/1943_2.png) [@mary\_h](https://fenicsproject.discourse.group/u/mary_h)\
**Post date:** [January 23, 2021, 9:37pm UTC](https://fenicsproject.discourse.group/t/derivative-function-has-no-attribute-subs/4956/6 "2021-01-23T21:37:05Z")

</div>

Hello,  
Regarding the code that you have provided, I was wondering what is the initial u vector that `assemble_system(jacobian, -residual, bcs)` uses to calculate the `b` and `res_initial` before entering the while loop?  
Also, did you change the `y` vector to the `u` vector on purpose?

Thanks a lot!

---

<div class="post-metadata">

**Author:** ![bhaveshshrimali](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bhaveshshrimali/32/124_2.png) [@bhaveshshrimali](https://fenicsproject.discourse.group/u/bhaveshshrimali)\
**Post date:** [January 24, 2021, 4:27pm UTC](https://fenicsproject.discourse.group/t/derivative-function-has-no-attribute-subs/4956/7 "2021-01-24T16:27:41Z")

</div>

It assembles the residual with the initial `u` (which in this case is simply zeros) after applying the dirichlet boundary conditions `bcs` and, in this case, the body force term (`B`). For this specific case, it is simply the contribution from the last term

I changed the formulation to the displacement field `u` since that is more common practice, but should work for the deformation field too. The example that I linked nicely illustrates the Newton-iteration loop embedded within a global step-loop.