# Results differ between dolfin and dolfinx

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)

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

def boundary_left_top(x):
return abs(x+H/2) < tol and (x)> tol

def boundary_left_bottom(x):
return abs(x+H/2) < tol and (x)< 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)

# Define strain-rate tensor
# Strain tensor
def epsilon(u):

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\

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 import mesh, fem, nls, io, la, log
import ufl
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.geo.synchronize()

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]
####################################################

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):

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!

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.

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

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
``````

``````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

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

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?