Difference in solution by UFL and tensorial notation for thermo-elasticity

1

Hi all,

I wish to solve transient coupled thermo-elasticity initial boundary value problem on circular domain with plane strain. I am going to complicate the model on top of this in later stages. For simplicity, I have provided a constant eigen-strain and solved for a single time increment. The variational formulation is written in both tensorial form and voigt form. Later on I need to use voigt form extensively. I am comparing the solution provided by both of them. They are a bit different (around 8%). This is not be significant for this simple case. But this error keeps on increasing for my complicated model. Below is my FEniCS code:

from dolfin import *
from fenics import *
import mshr
import os
set_log_level(ERROR)

dir_path = os.path.dirname(os.path.realpath(__file__))

######################################################################################################
R=10.00
dt = 0.001
E = 40
nu = 0.3
mu = E/(2.0*(1.0+nu)) # Lame's parameters
alpha = 0.1
lamb = E*nu/((1.0+nu)*(1.0-2.0*nu))
beta = E*alpha/((1.0-2.0*nu)*3.0)

############## MESH  ##################################
domain = mshr.Circle(Point(0,0),R) 
mesh = mshr.generate_mesh(domain,41)

####################################################
V = VectorElement("CG", mesh.ufl_cell(), 2) # Displacement u
Q = FiniteElement("CG", mesh.ufl_cell(), 2) # Temperature

a = MixedElement([V,Q])
W = FunctionSpace(mesh, a)
w = Function(W) 
(u, p) = split(w)
(v, q) = TestFunctions(W)
dw = TrialFunction(W)

class centre(SubDomain):		
    def inside(self,x, on_boundary):
     tol = 1e-14
     return (near(x[0],0.00) and near(x[1],0.00))

class circle_boundary(SubDomain):
    def inside(self,x, on_boundary):
      tol = 1e-1
      r=(x[0]**2+x[1]**2)**0.5
      return ( R-r < tol) and on_boundary



centre= centre()
circle_boundary = circle_boundary()



#DEFINING BOUNDARY CONDITIONS
bc1= DirichletBC(W.sub(0),Constant((0.00,0.00)),centre,method='pointwise') # Centre of circle is fixed
bc2=DirichletBC(W.sub(1),Constant((0.1)),circle_boundary) # on boundary temperature is applied
bc=[bc1,bc2]
d=2;

#DEFINING TERMS 
def epsilon_ufl(u):
			return as_vector([u[i].dx(i) for i in range(2)] +
                            [0.5*u[i].dx(j) + 0.5*u[j].dx(i) for (i,j) in [(0,1)]])
   
def sigma_ufl(u,p):
   cons = E/((1+nu)*(1-2*nu))
   #D_aniso = as_matrix([[cons*(1-nu),cons*nu,0.0],[cons*nu,cons*(1-nu),0.0],[0.0,0.0 ,cons*(1-2*nu)]])
   D_aniso = as_matrix([[2*mu+lamb, lamb,0.0 ],[ lamb,2*mu+lamb,0.0],[ 0.0,0.0 ,2*mu]])
   sigma_voigt = dot(D_aniso,epsilon_ufl(u))
   sigma1 = as_matrix([[sigma_voigt[0],sigma_voigt[2]],[sigma_voigt[2],sigma_voigt[1]]]) - 0.1*Identity(d)
   return sigma1



def epsilon(u):
 return 0.5*(nabla_grad(u)+nabla_grad(u).T)
 

def sigma(u,p):
 return lamb*tr(epsilon(u))*Identity(d)+2*mu*epsilon(u) -0.1*Identity(d)

  

n = 0

w0 =Function(W)
(u0, p_n) = split(w0)
p_n=project(Constant(0.99),W.sub(1).collapse())

##########################################################################
ptA = Point(10.0, 0.0)
n=0
F= inner(sigma(u,p),sym(grad(v)))*dx + dot(p,q)*dx + dot(dt,inner(grad(p), grad(q)))*dx - dot((p_n),q)*dx 
F_ufl = inner(sigma_ufl(u,p),sym(grad(v)))*dx + dot(p,q)*dx + dot(dt,inner(grad(p), grad(q)))*dx - dot((p_n),q)*dx 

for n in range(1):
    solve(F == 0, w, bc,solver_parameters={"newton_solver":{"relative_tolerance":1e-10,"absolute_tolerance":1e-10}}, form_compiler_parameters={"quadrature_degree":12})    
    print ('Tensorial notation', w(ptA))
    #solve(F_ufl == 0, w, bc,solver_parameters={"newton_solver":{"relative_tolerance":1e-10,"absolute_tolerance":1e-10}}, form_compiler_parameters={"quadrature_degree":12})    
    #print ('UFL', w(ptA))

When I run this code, the solution provided is

('Tensorial notation', array([ 0.01308619,  0.00024542,  0.1       ]))

And when I comment out

solve(F == 0, w, bc,solver_parameters={"newton_solver":{"relative_tolerance":1e-10,"absolute_tolerance":1e-10}}, form_compiler_parameters={"quadrature_degree":12})    
print ('Tensorial notation', w(ptA))

and uncomment

solve(F_ufl == 0, w, bc,solver_parameters={"newton_solver":{"relative_tolerance":1e-10,"absolute_tolerance":1e-10}}, form_compiler_parameters={"quadrature_degree":12})    
print ('UFL', w(ptA))

the solution is

('UFL', array([ 0.01306909,  0.00026225,  0.1       ]))

Please notice the difference, especially in second component. The difference is increasing in my complicated model. Am I missing anything? Any pointers will be valuable. Thanks in advance

Dolfin-version: 2017.2.0
OS: Ubuntu 16.04 LTS

Regards,
Peter

3

I slightly modified your script and ran it with docker:
docker run -ti --rm -v $(pwd):/home/fenics/shared/ -w /home/fenics/shared quay.io/fenicsproject/stable:2017.2.0
I cant reproduce your error:

from dolfin import *
from fenics import *
import mshr
import os
set_log_level(ERROR)

dir_path = os.path.dirname(os.path.realpath(__file__))

######################################################################################################
R=10.00
dt = 0.001
E = 40
nu = 0.3
mu = E/(2.0*(1.0+nu)) # Lame's parameters
alpha = 0.1
lamb = E*nu/((1.0+nu)*(1.0-2.0*nu))
beta = E*alpha/((1.0-2.0*nu)*3.0)

############## MESH  ##################################
domain = mshr.Circle(Point(0,0),R) 
mesh = mshr.generate_mesh(domain,41)

####################################################
V = VectorElement("CG", mesh.ufl_cell(), 2) # Displacement u
Q = FiniteElement("CG", mesh.ufl_cell(), 2) # Temperature

a = MixedElement([V,Q])
W = FunctionSpace(mesh, a)
w = Function(W) 
(u, p) = split(w)
(v, q) = TestFunctions(W)
dw = TrialFunction(W)

class centre(SubDomain):		
    def inside(self,x, on_boundary):
     tol = 1e-14
     return (near(x[0],0.00) and near(x[1],0.00))

class circle_boundary(SubDomain):
    def inside(self,x, on_boundary):
      tol = 1e-1
      r=(x[0]**2+x[1]**2)**0.5
      return ( R-r < tol) and on_boundary
centre= centre()
circle_boundary = circle_boundary()

#DEFINING BOUNDARY CONDITIONS
bc1= DirichletBC(W.sub(0),Constant((0.00,0.00)),centre,method='pointwise') # Centre of circle is fixed
bc2=DirichletBC(W.sub(1),Constant((0.1)),circle_boundary) # on boundary temperature is applied
bc=[bc1,bc2]
d=2;

#DEFINING TERMS 
def epsilon_ufl(u):
			return as_vector([u[i].dx(i) for i in range(2)] +
                            [0.5*u[i].dx(j) + 0.5*u[j].dx(i) for (i,j) in [(0,1)]])
   
def sigma_ufl(u,p):
   cons = E/((1+nu)*(1-2*nu))
   #D_aniso = as_matrix([[cons*(1-nu),cons*nu,0.0],[cons*nu,cons*(1-nu),0.0],[0.0,0.0 ,cons*(1-2*nu)]])
   D_aniso = as_matrix([[2*mu+lamb, lamb,0.0 ],[ lamb,2*mu+lamb,0.0],[ 0.0,0.0 ,2*mu]])
   sigma_voigt = dot(D_aniso,epsilon_ufl(u))
   sigma1 = as_matrix([[sigma_voigt[0],sigma_voigt[2]],[sigma_voigt[2],sigma_voigt[1]]]) - 0.1*Identity(d)
   return sigma1

def epsilon(u):
 return 0.5*(nabla_grad(u)+nabla_grad(u).T)
 
def sigma(u,p):
 return lamb*tr(epsilon(u))*Identity(d)+2*mu*epsilon(u) -0.1*Identity(d)

  

n = 0
w0 =Function(W)
(u0, p_n) = split(w0)
p_n=project(Constant(0.99),W.sub(1).collapse())

##########################################################################
ptA = Point(10.0, 0.0)
n=0
F= inner(sigma(u,p),sym(grad(v)))*dx + dot(p,q)*dx + dot(dt,inner(grad(p), grad(q)))*dx - dot((p_n),q)*dx 
F_ufl = inner(sigma_ufl(u,p),sym(grad(v)))*dx + dot(p,q)*dx + dot(dt,inner(grad(p), grad(q)))*dx - dot((p_n),q)*dx 

for n in range(1):
    solve(F == 0, w, bc,solver_parameters={"newton_solver":{"relative_tolerance":1e-10,"absolute_tolerance":1e-10}}, form_compiler_parameters={"quadrature_degree":12})
    w1 = w.copy(deepcopy=True)
    print ('Tensorial notation', w(ptA))
    solve(F_ufl == 0, w, bc,solver_parameters={"newton_solver":{"relative_tolerance":1e-10,"absolute_tolerance":1e-10}}, form_compiler_parameters={"quadrature_degree":12})    
    print ('UFL', w(ptA))
from numpy import allclose
print(allclose(w.vector().get_local(), w1.vector().get_local()))

Which produces the output:

fenics@ba1f4ce6c9c1:~/shared$ python3 source.py 
Tensorial notation [ 0.01196359 -0.00440094  0.1       ]
UFL [ 0.01196359 -0.00440094  0.1       ]
True
4

If I add

w.assign(Function(W))

between the nonlinear solves, I reproduce the error. It seems that the solution from the first problem satisfies the second one to within solver’s absolute tolerance, but, if you solve again from a zero initial guess, you get a slightly different solution. (If you comment-out set_log_level(ERROR), which I happened to do for compatibility with version 2018.1, you will see that the second solve finished in 0 iterations and 0 linear solver iterations, so w is unchanged.)

5

Nice catch, as the ufl forms are identical,
verified by doing:
allclose(assemble(F_ufl).get_local(), assemble(F).get_local())
I guess the difference must be in the Jacobian.

6

Thanks for the reply, @kamensky and @dokken. As @kamensky noted rightly, the error is not reproduced because the solution from the first problem satisfies the second one to within solver’s absolute tolerance. Still I am clueless about the difference. Please can you help me? Problem still remains unsolved.

Thanks in advance!
Peter