Updating error in newton solver

variational formulation
1

I want to solve a mixed function space problem, I have two function space P: scaler space and u: vector space which are defined as:

U_ = TestFunction(V)
(u_, P_) = split(U_)
dU = TrialFunction(V)
(du, dP) = split(dU) 
U = Function(V)
(u, P) = split(U)

For solving the problem I defined a zero initial value for P as:

P_init = Expression ( "0.0", degree = 0 )
P_old = interpolate ( P_init, V.sub(1).collapse())

Since I want to solve a time dependent problem, I need to update P_old in each step of Newton solver. The governing equations are defined as:

mech_form = inner(F1, grad(u_))*dx + inner(140*f,u_)*dx
phi_form = J*inner(((par_psi_phi) + 6.5*f), grad(P_))*dx + ( P - P_old )/dt * P_ * dx - J*g*P_*ds(1)
F = mech_form + phi_form
J = derivative(F, U, dU)

The Newton solver has written as:

U = Function(V)
(u1, P1) = split(U)
for (i, dti) in enumerate(np.array(t)):
    dt.assign(dti)
    m = np.hstack((m,dti))
    (u1, P1) = split(U)
    P_old.assign(P1)
    solver.solve()

but I got the following error: Expected a Function or linear combinations of Functions in the same FunctionSpace

I would be a

2

Hi,
It seems like a part of your post is missing. Please supply the rest, as well as sufficient code to be able to make a running example that can produce the Error message. Right now your form is independent of dU and dP

3

Hi Dokken, thank you for your reply, I have attached my code here.

from fenics import *
from dolfin import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
import pylab  
import random
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

L = 0.5
R = 0.5
N = 60 # mesh density
domain = Rectangle(Point(-70.0,-100.0), Point(70.0, 40.0)) - Rectangle(Point(-2.0,0.0), Point(2.0, 40.95))#-Circle(Point(0., 0.), R,100)
mesh = generate_mesh(domain, N)
plot ( mesh, title = 'mesh' )
mesh_points=mesh.coordinates()

Dphi = Constant(1.0)
E = 1000000
nu = 0.3
lmbda = Constant(E*nu/((1+nu)*(1-2*nu)))
mu = Constant(E/2/(1+nu))

d = 1 # interpolation degree
Vue = VectorElement('CG', mesh.ufl_cell(), d) # displacement finite element
Vpe = FiniteElement('CG', mesh.ufl_cell(), d) # temperature finite element
V = FunctionSpace(mesh, MixedElement([Vue, Vpe]))

def inner_b(x, on_boundary):
    return near(x[1], -0.0) and on_boundary
def inner_l(x, on_boundary):
    return near(x[0], -2.0) and on_boundary
def inner_r(x, on_boundary):
    return near(x[0], 2.0) and on_boundary
def bottom(x, on_boundary):
    return near(x[1], -100.0) and on_boundary
def left(x, on_boundary):
    return near(x[0], -70.0) and on_boundary
def right(x, on_boundary):
    return near(x[0], 70.0) and on_boundary
def top(x, on_boundary):
    return near(x[1], 40.0) and on_boundary
bc1 = DirichletBC(V.sub(1), Constant(0.), left)
bc2 = DirichletBC(V.sub(1), Constant(0.), right)
bc3 = DirichletBC(V.sub(1), Constant(0.), bottom)
bc7 = DirichletBC(V.sub(0).sub(0), Constant(0.), left)
bc8 = DirichletBC(V.sub(0).sub(0), Constant(0.), right)
bc9 = DirichletBC(V.sub(0).sub(0), Constant(0.), bottom)
bc10 = DirichletBC(V.sub(0).sub(1), Constant(0.), left)
bc11 = DirichletBC(V.sub(0).sub(1), Constant(0.), right)
bc12 = DirichletBC(V.sub(0).sub(1), Constant(0.), bottom)
bcs = [bc1,bc2,bc3,bc7,bc8,bc9,bc10,bc11,bc12]

# Defining multiple Neumann boundary conditions 
mf = MeshFunction("size_t", mesh, 1)
mf.set_all(0) # initialize the function to zero
class inner_b(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0) and on_boundary
inner_b = inner_b() # instantiate it
inner_b.mark(mf, 1)
ds = ds(subdomain_data = mf)

U_ = TestFunction(V)
(u_, P_) = split(U_)
dU = TrialFunction(V)
(du, dP) = split(dU) 
U = Function(V)
(u, P) = split(U)

d = u.geometric_dimension()
I = Identity(d)             # Identity tensor

F = I + grad(u)             # Deformation gradient
C = F.T * F                # Elastic right Cauchy-Green tensor
Finv = inv(F)
k = 0.005**2   #   = ro * D 
phi = 0.1
R = 0.08206*101325
T = 298

# Invariants of deformation tensors
Ic = tr(C)
J = det(F)

psi = (1-phi)*((mu/2)*(Ic - tr(I)) - mu*ln(J) + (lmbda/2)*(ln(J))**2)-phi*P*R*T*(1+1.5*ln(1.5*R*T)-ln(P*1/J))

F1 = (1-phi)* (mu * F - mu * Finv.T + lmbda*ln(J) * Finv.T) - phi*(P*R*T* Finv.T) # S

par_psi_phi = k*R*T*phi*grad(P/J) #phi*R*T*(-1.5*ln(1.5*R*T)+ln(P*1/1))  # P_p

#  Define time things.
dt = Constant(0.01)

P_init = Expression ( "0.0", degree = 0 )
P_old = project ( P_init, V.sub(1).collapse())

y_BC = Expression(("0.0", "0.0"),degree=0)
u = project(y_BC,V.sub(0).collapse())

f = Expression(("0.0", "0.0"),degree=0) #Constant((0,-10))
g = Expression(("5.0"),degree=0)

mech_form = inner(F1, grad(u_))*dx + inner(140*f,u_)*dx
phi_form = 1*J*inner(((par_psi_phi) + 6.5*f), grad(P_))*dx + ( P - P_old )/dt * P_ * dx - J*g*P_*ds(1)
F = mech_form + phi_form
J = derivative(F, U, dU)

problem = NonlinearVariationalProblem(F, U, bcs, J)
solver  = NonlinearVariationalSolver(problem)

prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-5
prm['newton_solver']['maximum_iterations'] = 25
prm['newton_solver']['relaxation_parameter'] = 1.0

m = []
Nincr = 10
t = np.logspace(0 ,2, Nincr+1)
Nx = 50 
x = np.linspace(R, L, Nx)
T_res = np.zeros((Nx, Nincr+1))
U = Function(V)
(u1, P1) = split(U)
for (i, dti) in enumerate(np.diff(t)):
    print ("Increment " + str(i+1))
    dt.assign(dti)
    m = np.hstack((m,dti))
    P_old.assign(P1)  # or P_old.assign(U.sub(1))
    solver.solve()