Results differ between dolfin and dolfinx

dolfinx
1

Hello,

I am trying to translate some code in fenics to dolfinx and I am getting slightly different results between the two cases. I am using dolfinx-0.6.0

I was hoping to get some help regarding this in case there is something wrong with my implementation.

Here is the fenics code

from __future__ import print_function
from fenics import *
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
from ufl import nabla_div
from ufl import sign
from mshr import *
from numpy import inf

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

# Elastic parameters
E  = 2.1E5
nu = 0.3
mu    = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))

# Geometry
L=1.2 #length
H=1.0 #horizontal dimension
gc= 2.7
ld=5.00E-2 # characteristic lengths
a0=0.25
h0=((0.063**2)-(0.063/2)**2)**0.5
pc=0.50-a0-h0 #precrack length
at=pc+a0+h0 # total initial crack length

# Generate mesh
#crack with sharp tip
domain = Polygon([dolfin.Point(-H/2,-L/2),dolfin.Point(H/2,-L/2),dolfin.Point(H/2,L/2),dolfin.Point(-H/2,L/2),dolfin.Point(-H/2,0.063/2),\
dolfin.Point(-H/2+a0,0.063/2),dolfin.Point(-H/2+a0+h0,0.0025),dolfin.Point(-H/2+a0+h0+pc,0.0),dolfin.Point(-H/2+a0+h0,-0.0025),\
dolfin.Point(-H/2+a0,-0.063/2),dolfin.Point(-H/2,-0.063/2)])

h=ld/3.5
n_mesh=round(1.414/(2*h)*((1.2**2+1)**0.5/(2**0.5))) #mesh parameter valid only for ASTM model
mesh = generate_mesh(domain, n_mesh)
nn = FacetNormal(mesh)


#loads and boundary conditions

# Define boundary phase field
tol = 1E-14 #tolerance boundary locations

def boundary_left_top(x):
	return abs(x[0]+H/2) < tol and (x[1])> tol
	
def boundary_left_bottom(x):
	return abs(x[0]+H/2) < tol and (x[1])< tol

markers = MeshFunction('size_t', mesh, 2, mesh.domains())
dx = Measure('dx', domain=mesh, subdomain_data=markers)

V=VectorFunctionSpace(mesh,'CG',2) #Function Space displacement function
D=FunctionSpace(mesh,'CG',2) #Function Space phase field function
u=Function(V)
d=Function(D)

dx_vl = dx(metadata={"quadrature_degree": 3}) #reduced integration for volumetric locking

# Define strain-rate tensor
# Strain tensor
def epsilon(u):
    return 0.5*(nabla_grad(u)+nabla_grad(u).T)

lagr_u_0=(1-d)*(1-d)*(mu*inner(epsilon(u),epsilon(u)))*dx+\
lmbda/2*(1-d)*(1-d)*((nabla_div(u))**2)*dx_vl\
+ gc*((d*d)/(2*ld)+0.5*ld*dot(grad(d),grad(d)))*dx

vv=TestFunction(V)
du=TrialFunction(V)
# Compute first variation of Pi (directional derivative about u in the direction of v)
F_u= derivative(lagr_u_0, u, vv)
# Compute Jacobian of F
J_u= derivative(F_u, u, du)
# Define boundary conditions
bcu_bottom  = DirichletBC(V, (0.0,0.0), boundary_left_bottom)
bcu_top=DirichletBC(V.sub(1), (4e-6), boundary_left_top)
bcu = [bcu_bottom, bcu_top];
problem = NonlinearVariationalProblem(F_u,u,bcu,J_u)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 5e-2
prm['newton_solver']['relative_tolerance'] = 5e-2
prm['newton_solver']['maximum_iterations'] = 200
prm['newton_solver']['error_on_nonconvergence'] = False
prm['newton_solver']['relaxation_parameter'] = 0.1
solver.solve()
print(f'\nfenics result\n2-norm of u = {np.linalg.norm(u.vector().get_local(),ord=2)}')

and then here is the equivalent dolfinx code.

import gmsh
from dolfinx.io.gmshio import model_to_mesh, read_from_msh

from dolfinx import mesh, fem, nls, io, la, log
import ufl
from ufl import tr,Identity,grad,dot,nabla_div,sym,outer,inner,dx,nabla_grad
from mpi4py import MPI
import numpy as np
import os
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType


## Start geometry and mesh ##
L=1.2 #length
H=1.0 #horizontal dim
a0 = 0.25
h0 = (0.063**2 - (0.063/2)**2)**0.5
pc = 0.5-a0-h0  # pre-crack length
at = pc+a0+h0   # total initial crack length
lc = 15e-3
top_marker = 1
bottom_marker = 2
right_marker = 3
left_up_marker = 4
left_down_marker = 5
plate_marker = 6

proc = MPI.COMM_WORLD.rank
gmsh.initialize()
if proc == 0:
    gmsh.model.add('fracture_geometry')
    
    gmsh.model.geo.addPoint(-H/2,-L/2,0,lc,1)
    gmsh.model.geo.addPoint(H/2,-L/2,0,lc,2)
    gmsh.model.geo.addPoint(H/2,L/2,0,lc,3)
    gmsh.model.geo.addPoint(-H/2,L/2,0,lc,4)
    gmsh.model.geo.addPoint(-H/2,0.063/2,0,lc,5)
    gmsh.model.geo.addPoint(-H/2+a0,0.063/2,0,lc,6)
    gmsh.model.geo.addPoint(-H/2+a0+h0,0.0025,0,lc,7)
    gmsh.model.geo.addPoint(-H/2+a0+h0+pc,0.0,0,lc,8)
    gmsh.model.geo.addPoint(-H/2+a0+h0,-0.0025,0,lc,9)
    gmsh.model.geo.addPoint(-H/2+a0,-0.063/2,0,lc,10)
    gmsh.model.geo.addPoint(-H/2,-0.063/2,0,lc,11)
    
    gmsh.model.geo.addLine(1,2,1)
    gmsh.model.geo.addLine(2,3,2)
    gmsh.model.geo.addLine(3,4,3)
    gmsh.model.geo.addLine(4,5,4)
    gmsh.model.geo.addLine(5,6,5)
    gmsh.model.geo.addLine(6,7,6)
    gmsh.model.geo.addLine(7,8,7)
    gmsh.model.geo.addLine(8,9,8)
    gmsh.model.geo.addLine(9,10,9)
    gmsh.model.geo.addLine(10,11,10)
    gmsh.model.geo.addLine(11,1,11)
    
    gmsh.model.geo.addCurveLoop([1,2,3,4,5,6,7,8,9,10,11],1)
    
    pl = gmsh.model.geo.addPlaneSurface([1],1)
    gmsh.model.geo.synchronize()
    
    gmsh.model.addPhysicalGroup(1,[3],top_marker)
    gmsh.model.addPhysicalGroup(1,[1],bottom_marker)
    gmsh.model.addPhysicalGroup(1,[2],right_marker)
    gmsh.model.addPhysicalGroup(1,[4],left_up_marker)
    gmsh.model.addPhysicalGroup(1,[11],left_down_marker)
    gmsh.model.addPhysicalGroup(2,[1],plate_marker)
    
    #gmsh.model.mesh.setRecombine(2, pl)    # quadrilateral mesh
    gmsh.model.mesh.generate(2)
    
    gmsh.write(f'fracture_geometry.msh')

fracture_mesh, cell_tags, facet_tags = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, 2)
gmsh.finalize()

nn = ufl.FacetNormal(fracture_mesh)
####################################################

## Define Function spaces ##
V = fem.VectorFunctionSpace(fracture_mesh,('CG',2)) #Function Space Picard iterations
D = fem.FunctionSpace(fracture_mesh,('CG',2)) #Function Space Picard iterations

# Elastic parameters
E = 2.1E5
nu = 0.3
mu = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))
gc = 1.0
ld = 2.5E-2 # length scale crack
R=0.0

## Define boundary conditions ##
left_up_facet = facet_tags.find(left_up_marker)
left_down_facet = facet_tags.find(left_down_marker)

left_down_dofs_u = fem.locate_dofs_topological(V.sub(1), fracture_mesh.topology.dim-1, left_down_facet)   #y-direction only
left_down_bc_u = fem.dirichletbc(ScalarType(0.0), left_down_dofs_u, V.sub(1))
left_up_dofs_u = fem.locate_dofs_topological(V.sub(1), fracture_mesh.topology.dim-1, left_up_facet) #y-direction only
left_up_bc_u = fem.dirichletbc(ScalarType(4e-6), left_up_dofs_u, V.sub(1))
bcu = [left_down_bc_u, left_up_bc_u]
####################################################

dx_vl = ufl.dx(metadata={"quadrature_degree": 3}) #reduced integration for volumetric locking
jit_options = {"timeout":60}
u = fem.Function(V)
u.name = "displacement"
d = fem.Function(D)
d.name = "phase_field"

# Solve for displacement #
# Define strain-rate tensor
def epsilon(u):
    return 0.5*(nabla_grad(u)+nabla_grad(u).T)

lagr_u = (1-d)*(1-d)*(mu*inner(epsilon(u),epsilon(u)))*dx + lmbda/2*(1-d)*(1-d)*((nabla_div(u))**2)*dx_vl + gc*((d*d)/(2*ld)+0.5*ld*dot(grad(d),grad(d)))*dx

vv=ufl.TestFunction(V)
du=ufl.TrialFunction(V)

F_u = ufl.derivative(lagr_u, u, vv)
J_u = ufl.derivative(F_u, u, du)

problem = fem.petsc.NonlinearProblem(F_u,u,bcu,J_u,jit_options=jit_options)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD,problem)
solver.atol = 5e-2
solver.rtol = 5e-2
solver.max_it = 200
solver.report = True
solver.relaxation_parameter = 0.1
solver.solve(u) 

print(f'\ndolfinx result\n2-norm of u = {u.x.norm()}')

I am looking at the 2-norm of the solution vector in both cases, and the results are different. Here it is for fenics and dolfinx respectively.

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 6.128e+00 (tol = 5.000e-02) r (rel) = 1.000e+00 (tol = 5.000e-02)
  Newton iteration 1: r (abs) = 5.515e+00 (tol = 5.000e-02) r (rel) = 9.000e-01 (tol = 5.000e-02)
  Newton iteration 2: r (abs) = 4.963e+00 (tol = 5.000e-02) r (rel) = 8.100e-01 (tol = 5.000e-02)
  Newton iteration 3: r (abs) = 4.467e+00 (tol = 5.000e-02) r (rel) = 7.290e-01 (tol = 5.000e-02)
  Newton iteration 4: r (abs) = 4.020e+00 (tol = 5.000e-02) r (rel) = 6.561e-01 (tol = 5.000e-02)
  Newton iteration 5: r (abs) = 3.618e+00 (tol = 5.000e-02) r (rel) = 5.905e-01 (tol = 5.000e-02)
  Newton iteration 6: r (abs) = 3.256e+00 (tol = 5.000e-02) r (rel) = 5.314e-01 (tol = 5.000e-02)
  Newton iteration 7: r (abs) = 2.931e+00 (tol = 5.000e-02) r (rel) = 4.783e-01 (tol = 5.000e-02)
  Newton iteration 8: r (abs) = 2.638e+00 (tol = 5.000e-02) r (rel) = 4.305e-01 (tol = 5.000e-02)
  Newton iteration 9: r (abs) = 2.374e+00 (tol = 5.000e-02) r (rel) = 3.874e-01 (tol = 5.000e-02)
  Newton iteration 10: r (abs) = 2.137e+00 (tol = 5.000e-02) r (rel) = 3.487e-01 (tol = 5.000e-02)
  Newton iteration 11: r (abs) = 1.923e+00 (tol = 5.000e-02) r (rel) = 3.138e-01 (tol = 5.000e-02)
  Newton iteration 12: r (abs) = 1.731e+00 (tol = 5.000e-02) r (rel) = 2.824e-01 (tol = 5.000e-02)
  Newton iteration 13: r (abs) = 1.558e+00 (tol = 5.000e-02) r (rel) = 2.542e-01 (tol = 5.000e-02)
  Newton iteration 14: r (abs) = 1.402e+00 (tol = 5.000e-02) r (rel) = 2.288e-01 (tol = 5.000e-02)
  Newton iteration 15: r (abs) = 1.262e+00 (tol = 5.000e-02) r (rel) = 2.059e-01 (tol = 5.000e-02)
  Newton iteration 16: r (abs) = 1.135e+00 (tol = 5.000e-02) r (rel) = 1.853e-01 (tol = 5.000e-02)
  Newton iteration 17: r (abs) = 1.022e+00 (tol = 5.000e-02) r (rel) = 1.668e-01 (tol = 5.000e-02)
  Newton iteration 18: r (abs) = 9.197e-01 (tol = 5.000e-02) r (rel) = 1.501e-01 (tol = 5.000e-02)
  Newton iteration 19: r (abs) = 8.277e-01 (tol = 5.000e-02) r (rel) = 1.351e-01 (tol = 5.000e-02)
  Newton iteration 20: r (abs) = 7.450e-01 (tol = 5.000e-02) r (rel) = 1.216e-01 (tol = 5.000e-02)
  Newton iteration 21: r (abs) = 6.705e-01 (tol = 5.000e-02) r (rel) = 1.094e-01 (tol = 5.000e-02)
  Newton iteration 22: r (abs) = 6.034e-01 (tol = 5.000e-02) r (rel) = 9.848e-02 (tol = 5.000e-02)
  Newton iteration 23: r (abs) = 5.431e-01 (tol = 5.000e-02) r (rel) = 8.863e-02 (tol = 5.000e-02)
  Newton iteration 24: r (abs) = 4.888e-01 (tol = 5.000e-02) r (rel) = 7.977e-02 (tol = 5.000e-02)
  Newton iteration 25: r (abs) = 4.399e-01 (tol = 5.000e-02) r (rel) = 7.179e-02 (tol = 5.000e-02)
  Newton iteration 26: r (abs) = 3.959e-01 (tol = 5.000e-02) r (rel) = 6.461e-02 (tol = 5.000e-02)
  Newton iteration 27: r (abs) = 3.563e-01 (tol = 5.000e-02) r (rel) = 5.815e-02 (tol = 5.000e-02)
  Newton iteration 28: r (abs) = 3.207e-01 (tol = 5.000e-02) r (rel) = 5.233e-02 (tol = 5.000e-02)
  Newton iteration 29: r (abs) = 2.886e-01 (tol = 5.000e-02) r (rel) = 4.710e-02 (tol = 5.000e-02)
  Newton solver finished in 29 iterations and 29 linear solver iterations.

fenics result
2-norm of u = 0.000283609098810126


INFO:    underlay of /etc/localtime required more than 50 (98) bind mounts
Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 10%] Meshing curve 2 (Line)
Info    : [ 20%] Meshing curve 3 (Line)
Info    : [ 30%] Meshing curve 4 (Line)
Info    : [ 40%] Meshing curve 5 (Line)
Info    : [ 50%] Meshing curve 6 (Line)
Info    : [ 60%] Meshing curve 7 (Line)
Info    : [ 70%] Meshing curve 8 (Line)
Info    : [ 80%] Meshing curve 9 (Line)
Info    : [ 90%] Meshing curve 10 (Line)
Info    : [100%] Meshing curve 11 (Line)
Info    : Done meshing 1D (Wall 0.00224078s, CPU 0s)
Info    : Meshing 2D...
Info    : Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.171377s, CPU 0.160911s)
Info    : 6766 nodes 13541 elements
Info    : Writing 'fracture_geometry.msh'...
Info    : Done writing 'fracture_geometry.msh'

dolfinx result
2-norm of u = 0.0003201253598044586
[WARNING] yaksa: 2 leaked handle pool objects

Thank you for any help I can get with this!

2

Your nonlinear convergence tolerances (atol and rtol) aren’t very tight, and as such your nonlinear solutions aren’t that well-converged. You’re comparing two 2-norms that are around 3e-4 while the tolerances of your newton solver are only 5e-2. I’d try running the same codes with tolerances that are (at the very least) significantly smaller than the norms of your solution vectors. Try 1e-6 or smaller to see whether the difference persists. Even just running the codes with a tolerance of 1e-4 (or similar) should already decrease the difference between both solution norms if there are no errors in your code.

3

Thank you for your input! I’ve re-run the codes with a tolerance (atol and rtol) of 1e-8, but the results are still very similar to what I initially reported.

fenics result
2-norm of u = 0.0002976277458932075

dolfinx result
2-norm of u = 0.000323976101263422

Thank you

4

I’m running a different version of Dolfinx so I’m not able to run your code. However, looking at your code snippets I see some differences in the parameter values used in both codes. In your Dolfin code:

gc= 2.7
ld=5.00E-2 # characteristic lengths

And in your Dolfinx code:

gc = 1.0
ld = 2.5E-2 # length scale crack

I’d also check whether both meshes are (near-)identical, if you haven’t done so already.

1 Like
5

Oh shoot, thank you for catching that. I made a couple of mistakes while copying stuff around. After I fixed that there was still some discrepancy in the norm values, which was closed when I used the exact same mesh file for both of them. I thought that the mesh was identical enough before, but apparently not!

1 Like
6

Hello, teja781,
I’m in the process of migrating my fenics phase field fracture codes to fenicsx, but I don’t have a clue about the phase field part, can you help me?