Hi all
I’m solving nonlinear poroelasticity problem. And same to the hyperelasticity demo, I use the minimization principle to generate the weak form F, then calculate the displacement field u and pressure field p by NonlinearVariationSolver for each time step. I read similar topics here and here, and follow what they are using. But I find for each time step, it always gives me the same result. I think it’s probably because the F isn’t updated even if I try to use “FunctionAssigner” to update the previous displacement field u_old. Here is my simplified code:
V = VectorElement(‘CG’,mesh.ufl_cell(),2)
Q = FiniteElement(‘CG’,mesh.ufl_cell(),1)
W = FunctionSpace(mesh,MixedElement([V,Q]))
w = Function(W, name=“Variables at current step”)
(u, p) = split(w)
w_old = Function(W, name=“Variables at previous step”)
(u_old, p_old) = split(w_old)
w_ = TestFunction(W)
(u_,p_,) = split(w_)
dw = TrialFunction(W)
(du,dp) = split(dw)#symmetric boundary conditions on left and bottom surfaces
bc1 = DirichletBC(W.sub(0).sub(0), Constant((0.0)), left)
bc2 = DirichletBC(W.sub(0).sub(1), Constant((0.0)), bottom)
bc3 = DirichletBC(W.sub(1),Constant((pb)),top)
bc4 = DirichletBC(W.sub(1),Constant((pb)),right)
bcs = [bc1, bc2, bc3, bc4]#initial conditions
uNot = Expression((‘lmbda0x[0]-x[0]', 'lmbda0x[1]-x[1]’), lmbda0 = lmbda0, degree = 2)
pNot = Constant(ln((pow(lmbda0,2)-1)/pow(lmbda0,2))+pow(lmbda0,-2)+xpow(lmbda0,-4)+Nv(1-pow(lmbda0,-2)))
print(ln((pow(lmbda0,2)-1)/pow(lmbda0,2))+pow(lmbda0,-2)+xpow(lmbda0,-4)+Nv(1-pow(lmbda0,-2)) )
#previous step variables
u_old = interpolate(uNot,W.sub(0).collapse())
p_old = interpolate(pNot,W.sub(1).collapse())#weak formulation
f = I + grad(u)
f0 = I + grad(u_old)
c = f.Tf
ic, j = tr(c), det(f)
j0 = det(f0)
m =(j- 1)dot(inv(f),(inv(f)).T)
free_energy=0.5(ic-3-2ln(j))+((j-1)ln((j-1)/j)+x(j-1)/j)/Nv
chemical_part=p*(j-j0)/Nv
dissipation_potential2=1/2/Nvinner(m,outer(grad(p),grad(p)))
incremental_potential=free_energydx -chemical_partdx-dtdissipation_potential2*dx
F = derivative(incremental_potential,w,w_)
dF = derivative(F, w, dw)
problem = NonlinearVariationalProblem(F,w,bcs,dF)
solver = NonlinearVariationalSolver(problem)#time loop
while ti < 1:
assign(w.sub(0),u_old)
assign(w.sub(1),p_old)
(no_of_iterations,converged)=solver.solve()
assigner = FunctionAssigner(W,W)
assigner.assign(w_old,w) #update w_old
(u_old, p_old) = w_old.split()
And if I put the following inside the time loop, it works. But it needs to calculate all the things for every time step, so the code will be much slower.
f = I + grad(u)
f0 = I + grad(u_old)
c = f.Tf
ic, j = tr(c), det(f)
j0 = det(f0)
m =(j- 1)dot(inv(f),(inv(f)).T)
free_energy=0.5(ic-3-2ln(j))+((j-1)ln((j-1)/j)+x(j-1)/j)/Nv
chemical_part=p*(j-j0)/Nv
dissipation_potential2=1/2/Nvinner(m,outer(grad(p),grad(p)))
incremental_potential=free_energydx -chemical_partdx-dtdissipation_potential2*dx
F = derivative(incremental_potential,w,w_)
dF = derivative(F, w, dw)
problem = NonlinearVariationalProblem(F,w,bcs,dF)
solver = NonlinearVariationalSolver(problem)
Any help and suggestion would be greatly appreciated.
Thank you!