Hello,
I am unable to reproduce the same value for the norm of a vector when I run the code in serial and parallel. Perhaps there is some flaw in my implementation? I appreciate any help I can get.
Here is the minimal working example. The code was run on dolfinx-0.6.0
import gmsh
from dolfinx.io.gmshio import model_to_mesh
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
## Define geometry and generate mesh ##
L=1.2 #length
H=1.0 #horizontal dim
a0 = 0.25
h0 = (0.063**2 - (0.063/2)**2)**0.5
lc = 7.5e-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.0,0,lc,7)
gmsh.model.geo.addPoint(-H/2+a0,-0.063/2,0,lc,8)
gmsh.model.geo.addPoint(-H/2,-0.063/2,0,lc,9)
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,1,9)
gmsh.model.geo.addCurveLoop([1,2,3,4,5,6,7,8,9],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,[9],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('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(0.01), 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.solve(u)
print(f'\nRank {MPI.COMM_WORLD.rank} : 2-norm of u = {u.x.norm()}')
And the results I got when I ran the code in serial is
Rank 0 : 2-norm of u = 3.5633605875906254
vs. the result when run on 2 cpus is
Rank 0 : 2-norm of u = 2.3389230164029495
Rank 1 : 2-norm of u = 2.3389230164029495
Thank you!
