time-dependent Nonlinear Mixed Method fails to update solution after first solve

Question

I am attempting to solve a time-dependent nonlinear mixed formulation of a system of equations. For some reason, after the first time step, all subsequent solutions are identical, even though it should be evolving in time. The mixed space is a scalar field (phi) and a vector field (u). The system has an initial condition for phi ('phi_k') and should evolve phi in time. The vector field u is fully determined at each time step, given the previous iteration for phi.
My scheme is essentially:

Vv = VectorFunctionSpace(mesh, "Lagrange", 2)
Vf = FunctionSpace(mesh, "Lagrange", 1)
VV = Vf*Vv

# Dirichlet BC applied to one subspace
bc = DirichletBC(VV.sub(1), u_0, boundary) 

# Define functions
(phi, u) = TrialFunctions(VV) #incremental phase, displacement
(v,tau) = TestFunctions(VV)  
w  = Function(VV)           # Solution from previous iteration
F = (all terms of the system of equations added)

while t <= T:
    R  = action(F, w) 
    DR  = derivative(R,  w)   # Gateaux derivative 
    problem = NonlinearVariationalProblem(R, w, bc, DR)
    solver  = NonlinearVariationalSolver(problem)
    solver.solve()
    (phi_s,u_s) = split(w)      #the solutions

    phi_k = project(phi_s, Vf)
    t += dt

A minimal working model is below:

from dolfin import *
import numpy as np

##############################
# Functions ##################
##############################
def boundary(x, on_boundary):
    # define the Dirichlet boundary
    return on_boundary

##############################
# System Params ##############
##############################
# Material/System Parameters
kappa = 1.0
chi = 1.0
E , nu = 10000.0, 0.40 #J/m^3, unitless
mu, lamb = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
Ec = Constant(0.05*mu)

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

# Timing
t = 0.0
dt = 0.01
T = 10*dt


##############################
# Mesh Params ################
##############################
# Create mesh and define function space
nx = 100
mesh = UnitSquareMesh(nx, nx)
Vv = VectorFunctionSpace(mesh, "Lagrange", 2)
Vf = FunctionSpace(mesh, "Lagrange", 1)
VV = Vf*Vv

# Dirichlet BC
xc, yc = 0.5, 0.5
U = 0.05 # =w/R
u_0 = Expression(('U*sqrt( pow(x[0]-xc,2)+pow(x[1]-yc, 2) )*cos(atan2(x[1]-yc,x[0]-xc))' ,
                  'U*sqrt( pow(x[0]-xc,2)+pow(x[1]-yc, 2) )*sin(atan2(x[1]-yc,x[0]-xc))'), U=U, xc=xc, yc=yc)
#apply to second space (Vv)
bc = DirichletBC(VV.sub(1), u_0, boundary) 


##############################
# Define variational problem #
##############################
# Define functions
(phi, u) = TrialFunctions(VV) #incremental phase, displ
(v,tau) = TestFunctions(VV) 
w  = Function(VV)           # Solution from previous iteration

# Initial condition for phase 
class initialPhase(Expression):
  def eval(self, value, x):
    xslit = 0.5
    yslit = 0.5
    W,a = phase_size()
    if (abs(x[0]-xslit) <= W) and (abs(x[1]-yslit) <= a):
        value[0] = abs(x[0]-xslit)/W*0.95 +0.05 # linear/abs
    elif (abs(x[0]-xslit) <= W) and (abs(x[1]-yslit) >= a) and ((abs(x[1]-yslit)-a)**2+(x[0]-xslit)**2 <= W**2):
        value[0] = sqrt( (abs(x[0]-xslit)/W)**2 + ((abs(x[1]-yslit)-a) /W)**2)*0.95 +0.05 #linear/abs
    else:
        value[0] = 1.0

def phase_size():
    W = 1.5 * 0.05 
    a = 5 * 0.05
    return W,a

phik0 = initialPhase()
# Projecting phi_k -- the initial condition for phi 
phi_k = interpolate(phik0, Vf)


# functions for variational problem
d = u.geometric_dimension()
epsilon = sym(nabla_grad(u))
sigma = 2*mu*(epsilon) + lamb*tr(epsilon)*Identity(d)
Es = 0.5*lamb*tr(epsilon)*tr(epsilon) + mu*inner(epsilon.T,epsilon)
g = 4*phi_k**3 - 3*phi_k**4
gprime = 12*(phi_k**2-phi_k**3)

# Force balance (a(u,v) = L(v) = Integrate[f*v, dx] )
au = -g*inner(sigma, sym(nabla_grad(tau)))*dx + gprime*inner(sigma, outer(grad(phi_k),tau) )*dx #using symmetry properties of sigma
Lu = inner(Constant((0.0,0.0)),tau)*dx

# dphi/dt
print('Setting up Linear problem for Phi-- Implicit Euler')
aphi = phi*v*dx + dt*chi*kappa* dot(nabla_grad(phi), nabla_grad(v)) *dx + dt*chi*( gprime *Es )*v*dx
fphi = phi_k  + dt*chi*( gprime *(Ec) )
Lphi = fphi*v*dx


# Mixed Method
amix = au + aphi
Lmix = Lu + Lphi

F = amix - Lmix

ind = 0

while t <=T:
    print 'time = ', t
    R  = action(F, w) 
    DR  = derivative(R,  w)   # Gateaux derivative 
    problem = NonlinearVariationalProblem(R, w, bc, DR)
    solver  = NonlinearVariationalSolver(problem)
    solver.solve()
    (phi_s,u_s) = split(w)

    phi_k = project(phi_s, Vf)
    t += dt
    ind +=1

# Plot the displacement solution
plot(u_s, mode = "displacement", interactive = True)

# Plot the phase solution
plot(phi_s, title="phase", interactive = True)

Thanks to anyone who has ideas here!

Comments

Just a guess, but it seems that you are not updating F inside the time-loop. You define it once at the start, but then it does not get updated anymore inside the loop. What you are doing may work, if the forms are working with pointers on phi_k, unfortunately I am not sure about how this is implemented. Anyhow, I always define my variational probelem inside the loop. Hope this helps.

---

Unfortunately, that does not work. Thanks for the suggestion, though.

Answer 1

Magne Nordaas has the right idea, but we need to take the subspace of the mixed space for the degrees of freedom to match. Thanks Magne!

assigner = FunctionAssigner(Vf, VV.sub(0))
assigner.assign(phi_k, phi_s)

Answer 2

This seems to also work, it seems:

phi_s, u_s = w.split()
phi_sV = interpolate(phi_s,Vf)
phi_k.vector()[:] = phi_sV.vector()

Answer 3

Use FunctionAssigner to update phi_k inside the time loop.

assigner = FunctionAssigner(Vf, VV)
assigner.assign(phi_k, phi_s)