Solution is not updated after the first time step in time-dependent nonlinear problem

variational formulation
1

Dear FEniCS Community,

I am a new user and this is my first time posting a question here: I apologize in advance for any mistakes or if my questions are too basic.
I am trying to use FEniCS to solve a time-dependent hyperelasticity problem, using a mixed formulation in which the unknowns are the displacement (a vector field u) and the pressure (a scalar field p). Starting from the hyperelasticity demo, I attempted to write my own code, but I have an issue: the program seems to correctly solve the problem at the first time step, but then nothing else happens, and the solution is not updated in the following iterations of the time loop.

Here is a simplified version of the code:

import dolfin
from dolfin import * 

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

from ufl import grad as ufl_grad
def Grad(v):
        return ufl_grad(v)

mesh = UnitCubeMesh(6, 6, 6) 

#Elements and Function Spaces
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)   #displacement u
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)   #pressure p
TH = MixedElement([P2, P1])                          #u + p
V = FunctionSpace(mesh, TH)             #(u,p) 

# Mark boundary subdomains 
upp = CompiledSubDomain("near(x[2], side) && on_boundary", side = 1.0) 
down = CompiledSubDomain("near(x[2], side) && on_boundary", side = 0.0) 

#Define parameters
u0 = Constant((0.0, 0.0, 0.0))
pp = Constant(0.0)
T = 1e-08
num_steps = 100 
dt = T / num_steps
TT  = Constant((0.0, 0.0, -0.1))  # Traction force on the boundary - current configuration
k = Constant(5e-9)
K0 = as_matrix([  [k, 0, 0],  [0, k, 0],  [0, 0, k]  ])       

#Boundary Conditions: zero displacement at z = 0, zero pressure at z = 1
bcd_z = DirichletBC(V.sub(0), u0, down)
bcp = DirichletBC(V.sub(1), pp, upp) 
bcs = [bcd_z, bcp] 

# Define functions
dup = TrialFunction(V)             # Incremental displacement and pressure
du, dp = split(dup) 

u_, p_ = TestFunctions(V) 

up  = Function(V)         # Unknown Displacement and Pressure (n+1) 
u, p = split(up) 

up_prev = Function(V) 			   # Previous iteration 
u_prev, p_prev = split(up_prev) 

up_init = Expression( ("0.0", "0.0", "0.0", "0.0"), degree=1)
up_prev = interpolate(up_init, V)   #Initial condition for u and p

# Kinematics     
d = len(u) 
I = Identity(d) 
F = I + Grad(u)  #Deformation gradient 

#Jacobian of F
def Js(u): 
        return det(F) 

#First Piola-Kirchhoff stress tensor - Mooney Rivlin 
def P(u, p): 
	C = (F.T)*F 
	Ic = variable(tr(C))
	IIc = variable(0.5*(tr(C)**2 - tr(C*C)))
	#Mooney-Rivlin strain energy 
	mu1 = Constant(0.4)
	mu2 = Constant(0.25)
	psi = (mu1/2)*(Ic - 3) + (mu2/2)*(IIc - 3)
	
	gamma1 = diff(psi, Ic) + Ic*diff(psi, IIc)
	gamma2 = -diff(psi, IIc)
	
	P_s = 2*F*(gamma1*I + gamma2*C) 
	 
	return P_s - Js(u)*p*inv(F).T 

#Permeability tensor
def K_star(K0, u): 
	return Js(u)*inv(F)*K0*inv(F).T 

#Traction in material configuration
def T_star(TT, u): 
	return Js(u)*inv(F).T*TT

dolfin.parameters["form_compiler"]["quadrature_degree"] = 6

#Loop for time stepping
t = 0 
n = 1 
while (t <= T): 
	t += dt 
	print("Iteration", n, "-th", "Time", t)	
	
	#Variational form 
	L = Js(u)*p_*dx + dt*inner(Grad(p_), K_star(K0, u)*Grad(p))*dx - inner(P(u, p), Grad(u_))*dx - Js(u_prev)*p_*dx + inner(T_star(TT, u), u_)*ds 
	j = derivative(L, up, dup)
	
	problem = NonlinearVariationalProblem(L, up, bcs, J=j) 
	solver = NonlinearVariationalSolver(problem) 
	solver.solve() 
	
	assigner = FunctionAssigner(V, V) 
	assigner.assign(up_prev, up) 
	 
	n = n+1

It looks as if the nonlinear variational form L (or its jacobian) is not updated properly after the first time step: the first iteration is fine and the solution behaves as expected, but then it remains constant. I tried to make use of “FunctionAssigner” instead of the simple “assign” command as suggested in a similar topic I found here, but unfortunately it did not solve my problem. Defining everything inside the time loop does not seem to fix the issue as well.

Any help or suggestion would be greatly appreciated. Thank you all in advance!

2

What do you expect to happen with the solution for a constant exterior force after the equilibrium has been found (for a rate-independent material model, that is)?

If you just write
TT = Expression(("0.0","0.0","A*t"), t=0, A=1./T, element = P2)
and in your loop after t=t+dt write
TT.t = t
the solution will be different in the next timestep since the exterior force will increase with time.

As a side note: You can define forms and solver outside of the time loop, so the lines

L = Js(u)*p_*dx + dt*inner(Grad(p_), K_star(K0, u)*Grad(p))*dx - inner(P(u, p), Grad(u_))*dx - Js(u_prev)*p_*dx + inner(T_star(TT, u), u_)*ds 
j = derivative(L, up, dup)
	
problem = NonlinearVariationalProblem(L, up, bcs, J=j) 
solver = NonlinearVariationalSolver(problem)

can all be placed before the loop.

3

Thank you for the answer and for the hint: now the solution is updated. I was trying to find the mistake in the code and I lost sight of the physical behaviour of the system!