Computing residual directly from energy

[bhaveshshrimali]

Consider for instance a usual hyperelasticity problem (#398) which is basically finding the minimizer of a scalar potential energy (Pi): $\arg\min_u \Pi (u)$

It would be useful to know constructs to be able to directly calculate the residual from the energy using automatic differentiation (The residual to jacobian part is already taken care of by Gridap). This is particularly useful in more complicated constitutive relations, where calculating the residual by hand would be time-consuming and prone to algebraic errors.

MWE

using Gridap
using LineSearches: BackTracking
E = 10.
ν = 0.3
const μ = E/2/(1+ν)
const λ = 2*μ*ν/(1-2*ν)
F(∇u) = one(∇u) + ∇u
J(F) = det(F)
B(x) = VectorValue(0.0, -0.5, 0.0)
T(x) = VectorValue(-0.1, 0.0, 0.0)
ψ(F) = μ * (0.5 * (F ⊙ F - 3.) - log(J(F))) + λ/2. * (J(F) - 1.)^2
S(F) = ∇(ψ)(F)
energy(u) = ∫( ψ∘F∘∇(u) )*dΩ
res(u, v) = ∫( (S∘F∘∇(u)) ⊙ ∇(v) - B ⊙ v)*dΩ 
domain = (0.0, 1.0, 0.0, 1.0, 0.0, 1.0)
partition = (5, 5, 5)
model = CartesianDiscreteModel(domain, partition)

labels = get_face_labeling(model)
add_tag_from_tags!(labels, "rightFace", [2, 4, 6, 8, 14, 16, 18, 20, 26])
add_tag_from_tags!(labels, "leftFace", [1, 3, 5, 7, 13, 15, 17, 19, 25])

reffe = ReferenceFE(lagrangian, VectorValue{3, Float64}, 1) 
V = TestFESpace(model, reffe, conformity=:H1, dirichlet_tags = ["leftFace", "rightFace"]) 

Ω = Triangulation(model)
degree = 2
dΩ = Measure(Ω, degree)

nls = NLSolver(
show_trace=true,
method=:newton,
linesearch=BackTracking()
)

g0(x) = VectorValue(0.0, 0.0, 0.0)
g1(x) = VectorValue(
    0.0,
    (0.5 + (x[2]-0.5) * cos(π/3) - (x[3] - 0.5) * sin(π/3) - x[2])/2,
    (0.5 + (x[2]-0.5) * sin(π/3) + (x[3] - 0.5) * cos(π/3) - x[3])/2
)

U = TrialFESpace(V, [g0, g1])

#FE problem
op = FEOperator(res, U, V)
solver = FESolver(nls)

x0 = zeros(Float64, num_free_dofs(V))
uh = FEFunction(U, x0)
uh, = solve!(uh, solver, op)

Problem:

How to calculate the gradient of the energy? If I change the definition of residual to

function S(∇u)
    μ * (F(∇u) - inv(F(∇u))') + λ * (J(F(∇u))) * (J(F(∇u)) - 1.) * inv(F(∇u))'
end

function res(u, v)
    ∫( (S∘∇(u)) ⊙ ∇(v) - B ⊙ v )*dΩ - ∫( T ⊙ v )*dΓ_t
end

everything works fine as it should. Also if I calculate the energy from the solution uh, namely

julia> sum(energy(uh))
0.11957195697992658

I get the correct result. It is just that the calculation of the residual from the energy using AD that is proving to be tricky.

Sorry for a bit long post, but I hope that the issue is clear :)

[fverdugo]

Hi @bhaveshshrimali

S(F) = ∇(ψ)(F) This does not work since the AD machinery is for functions taking Point i.e. VectorValue as argument, not 2nd order tensors.

To auto generate a residual and a jacobian you need to do something like:

energy(u) = ∫( ψ∘F∘∇(u) )*dΩ # You also will need to add the external loads here
res(u,v) = gradient(energy,u)
jac(u,du,v) = hessian(energy,u) # I have realized that you have to explicitly build the jacobian like this
op = FEOperator(res, jac, U, V)

I have done some tests and it seems that the compilation of the jacobian is EXTREMELY slow, several minutes. The compilation of the residual is much more acceptable.

[fverdugo]

BTW, μ and λ are not defined in the MWE. Can you update the code with some values?

[bhaveshshrimali]

res(u,v) = gradient(energy,u)
jac(u,du,v) = hessian(energy,u) 

Thanks a lot @fverdugo ! Indeed this seems to work and not throw the error anymore. Is building the jacobian jac explicitly needed? Or can I just leave it to be computed under the hood. And is that different from this? I guess not, right?

A side question: how are you computing the individual time for the jacobian compilation? TimerOutputs? And is this (the slowness) an issue with ForwardDiff?

[bhaveshshrimali]

BTW, μ and λ are not defined in the MWE. Can you update the code with some values?

Done.

[fverdugo]

A side question: how are you computing the individual time for the jacobian compilation?

I explicitly call jac.

u = zero(U)
v = get_cell_shapefuns(V)
du = get_cell_shapefuns_trial(U)
j = jac(u,du,v)

It takes ages.

r = res(u,v)

is much faster

And is this (the slowness) an issue with ForwardDiff?

Possibly, but I don't know.