I have a highly nonlinear problem where the rhs depends on the solution, K(u)u=F(u). I am solving this by calling NonlinearVariationalSolver in a loop that recomputes and increments the rhs. For a MWE, I attach a slightly modified version of the Hyperelasticity demo that basically does the same thing. My problem is that, with this approach, each load step iteration is differentiated in the gradient computations. However, I just need the final solution of the nonlinear equation and not the series of nonlinear equations actually solved (apart from terrible performance, this is making it hard for me to debug other issues). In other words, I want solving a nonlinear equation in n amount of load steps to result in only 1 adjoint equation. From my understanding, I will probably have to overload a doflin-adjoint solve function and/or the SolveBlock so that all the load steps are treated as the outcome of one solve, but I am not sure where to start. Thanks for any suggestions in advance!

```
import matplotlib.pyplot as plt
from dolfin import *
from dolfin_adjoint import *
import numpy as np
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)
c = Constant((0.0, 0.0, 0.0))
r = Expression(("scale*0.0",
"scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
"scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3, degree=2)
bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr]
# Define functions
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
# Kinematics
I = Identity(3) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# Elasticity parameters
E, nu = Constant(10.0, name="E"), Constant(0.3)
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
# Total potential energy
Pi = psi*dx
T = Constant((0.5, 0.0, 0.0), name="traction") # Traction force on the boundary
def update_traction(i): # function to update load vector during solve can be arbitrarly complex
return Constant(i)*T
tape = get_working_tape()
file = File("displacement.pvd"); file << (u,-1)
for i in np.linspace(0, 1.0, 5): # 5 load steps
with tape.name_scope("Timestep"):
inc_load = update_traction(i)
Pi_f = dot(inc_load, u)*ds
F = derivative(Pi - Pi_f, u, v)
J = derivative(F, u, du)
problem = NonlinearVariationalProblem(F, u, bcs=bcs, J=J, form_compiler_parameters={"optimize": True})
solver = NonlinearVariationalSolver(problem)
solver.parameters.update({"nonlinear_solver":"snes"})
solver.solve()
file << (u,i)
tape.visualise()
J = assemble(psi*dx)
Jhat = ReducedFunctional(J, Control(mu))
h = Constant(0.01)
assert(taylor_test(Jhat, mu, h) > 1.9)
assert(taylor_test(Jhat, mu, h, dJdm=0) > 0.9)
```