Parallel implementation of adaptive mesh refinemet

Hello everyone,

I’m working with FEniCSx 0.8.0 to implement a brittle fracture phase field model using adaptive mesh refinement. The simulation behaves as expected when run in serial, but produces incorrect results when executed in parallel . I’ve been trying to debug the issue, but I’m not sure where to start, and I haven’t had much luck finding similar cases or solutions.

Below, I’ve attached my code. Sry for the lengthy code — it’s not a mwe since I’m importing the mesh from Gmsh. I’ve also attached the geo file used to generate the mesh.

import dolfinx.fem.petsc
import dolfinx.nls.petsc
from dolfinx import fem, io, log, mesh
import basix
import ufl
import numpy as np
from scipy.spatial import cKDTree
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
from mpi4py import MPI
              
MPI.COMM_WORLD.Barrier()
start_time = MPI.Wtime()    

meshfilename = "senp.msh"
if meshfilename.endswith(".msh"):
    dmesh, ct, ft = io.gmshio.read_from_msh(meshfilename, MPI.COMM_WORLD, 0, gdim=2)

P = basix.ufl.element("CG", dmesh.basix_cell(), 1, shape=(dmesh.geometry.dim,))
S = basix.ufl.element("CG", dmesh.basix_cell(), 1)
R = basix.ufl.element("DG", dmesh.basix_cell(), 0)
V = fem.functionspace(dmesh,P)
VW = fem.functionspace(dmesh, S)
W = fem.functionspace(dmesh,R)
du = ufl.TestFunction(V)
dd = ufl.TestFunction(VW)
HH = fem.Function(W)
Hold = fem.Function(W) # history field at the previous time step
u = ufl.TrialFunction(V)
d = ufl.TrialFunction(VW)
uold = fem.Function(V)
unew = fem.Function(V)
dold = fem.Function(VW)
dnew = fem.Function(VW)

E =  210e6
nu = 0.3
lmda = (E*nu) / ((1+nu)*(1-2*nu))
mu = E / (2*(1+nu))
K = (2.0/3.0)*mu + lmda
G_cr = 2.7e3
a=1 #mm
l = 0.01#linesearch
total_u = 0.02
delta_u = 0.5e-3  #for 0.3e-3
k = 1e-5

def top(x):
    return np.isclose(x[1],1)
def bot(x):
    return np.isclose(x[1],0)
def crack(x):
    return np.logical_and(abs(x[1]-0.5) <= 0.001, x[0] <= 0.5)

V01, _ = V.sub(0).collapse()
V02, _ = V.sub(1).collapse()
bot_facets = mesh.locate_entities_boundary(dmesh, 1, bot)
top_facets = mesh.locate_entities_boundary(dmesh, 1, top)
crack_facets = mesh.locate_entities_boundary(dmesh, 1, crack)
bot_dofs = fem.locate_dofs_topological(V, 1, bot_facets)
top_dofs = fem.locate_dofs_topological((V.sub(0),V01), 1, top_facets)
top_dofs_y = fem.locate_dofs_topological((V.sub(1),V02), 1, top_facets) 
crack_dofs = fem.locate_dofs_topological(VW, 1, crack_facets)
u_bot = fem.Function(V)
u_botval = fem.Constant(dmesh,ScalarType([0.0,0.0]))
u_botexpr = fem.Expression(u_botval,V.element.interpolation_points())
u_bot.interpolate(u_botexpr)
u_top = fem.Function(V01)
u_topval = fem.Constant(dmesh,ScalarType(1))    
u_topexpr = fem.Expression(u_topval,V01.element.interpolation_points())
u_top.interpolate(u_topexpr)
u_topy = fem.Function(V02)
u_topyval = fem.Constant(dmesh,ScalarType(0))    
u_topyexpr = fem.Expression(u_topyval,V02.element.interpolation_points())
u_topy.interpolate(u_topyexpr)
d_crack = fem.Function(VW)
d_crack.x.array[:] = 1 
dc = fem.Constant(dmesh, ScalarType(1))
bc_ubot = fem.dirichletbc(u_botval,bot_dofs,V)
bc_utop = fem.dirichletbc(u_top, top_dofs, V.sub(0))
bc_utopy = fem.dirichletbc(u_topy, top_dofs_y, V.sub(1))
bc_dcrack = fem.dirichletbc(dc, crack_dofs, VW)
t_bar = fem.Constant(dmesh,np.array([0.0,0.0],dtype=ScalarType))
b = fem.Constant(dmesh,np.array([0.0,0.0],dtype=ScalarType))
n = ufl.FacetNormal(dmesh)
ds_top = ufl.Measure(integral_type='ds', domain=dmesh, subdomain_data=ft, subdomain_id=2)

def invariants_principal(A):
    i1 = ufl.tr(A)
    i2 = (ufl.tr(A)**2 - ufl.tr(A * A)) / 2
    i3 = ufl.det(A)
    return i1, i2, i3
   
def eigenstate2(A):
    if ufl.shape(A) != (2, 2):
        raise RuntimeError(f"Tensor A of shape {ufl.shape(A)} != (2, 2) is not supported!")
    eps = 3.0e-16
    A = ufl.variable(A)
    I1, _, _ = invariants_principal(A)
    Δ = (A[0, 0] - A[1, 1])**2 + 4 * A[0, 1] * A[1, 0]  # = I1**2 - 4 * I2
    Δ += eps**2
    λ = (I1 - ufl.sqrt(Δ)) / 2, (I1 + ufl.sqrt(Δ)) / 2
    E = [ufl.diff(λk, A).T for λk in λ]
    return λ, E

def bracket_operator_plus(x):
    return 0.5*(x+abs(x))
def bracket_operator_minus(x):
    return 0.5*(x-abs(x))
def degradation_fn(x_d):
    return (1.0-x_d)**2
def spectral_positive_part(T):
    if ufl.shape(T) == (2, 2):
        eig, eig_vec = eigenstate2(T)
        T_p = bracket_operator_plus(eig[0]) * eig_vec[0]
        T_p += bracket_operator_plus(eig[1]) * eig_vec[1]
    return T_p


def spectral_negative_part(T):
    if ufl.shape(T) == (2, 2):
        eig, eig_vec = eigenstate2(T)
        T_p = bracket_operator_minus(eig[0]) * eig_vec[0]
        T_p += bracket_operator_minus(eig[1]) * eig_vec[1]
    return T_p
def epsilon(x_u):
    return ufl.sym(ufl.grad(x_u))
def sigma0(x_u):
    '''Miehe split'''
    eps = epsilon(x_u)
    eps_P = spectral_positive_part(eps)
    eps_N = spectral_negative_part(eps)
    return (lmda*bracket_operator_plus(ufl.tr(eps))*ufl.Identity(len(x_u)) + 2 * mu * eps_P) +\
        (lmda*bracket_operator_minus(ufl.tr(eps))*ufl.Identity(len(x_u)) + 2 * mu * eps_N)
def sigma(x_u, x_d):
    return (degradation_fn(x_d)+ k)*sigma0(x_u)
def sigmah(x_u):
    return 2.0*mu*epsilon(x_u) + lmda*ufl.tr(epsilon(x_u))*ufl.Identity(len(x_u))
def elastic_en(x_u, x_d):
    return 0.5 * ufl.inner(sigma(x_u, x_d), epsilon(x_u))*ufl.dx - ufl.dot(b,u)*ufl.dx - ufl.dot(t_bar,u)*ds_top
def frac_en(x_u, x_d):
    return G_cr * (x_d**2 / (2*l) + 0.5*l * ufl.dot(ufl.grad(x_d), ufl.grad(x_d)))*ufl.dx
def tot_en(x_u, x_d):
    total_energy= fem.assemble_scalar(fem.form((elastic_en(x_u, x_d) + frac_en(x_u, x_d)),jit_options={"timeout":50}))
    total_energy=dmesh.comm.allreduce(total_energy,op=MPI.SUM)
    return total_energy
def psi_e(x_u):
    eps = epsilon(x_u)
    eps_P = spectral_positive_part(eps)
    return 0.5 * lmda*bracket_operator_plus(ufl.tr(eps))**2 + mu * ufl.inner(eps_P, eps_P)
def H(x_uold,x_u,x_Hold):
    return ufl.conditional(ufl.operators.gt(psi_e(x_u),x_Hold),psi_e(x_u),x_Hold)

def obtain_midpoints(mesh):
    element_list = np.arange(mesh.topology.index_map(mesh.topology.dim).size_local, dtype=np.int32)
    mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
    midpoints = dolfinx.mesh.compute_midpoints(mesh, mesh.topology.dim, element_list).T
    return mesh , element_list, midpoints

def element_size_list(mesh):
    Size = dolfinx.fem.functionspace(mesh, ("DG", 0))  # Discontinuous Galerkin space (one value per element)
    h = ufl.CellDiameter(mesh)  # UFL symbolic representation of cell diameter
    h_func = dolfinx.fem.Function(Size)
    h_expr = dolfinx.fem.Expression(h, Size.element.interpolation_points())
    h_func.interpolate(h_expr)
    element_sizes = h_func.x.array
    return element_sizes

def filter_elements_by_size(mesh, element_ids, mesh_size):
    Size = dolfinx.fem.functionspace(mesh, ("DG", 0))  # Discontinuous Galerkin space (one value per element)
    h = ufl.CellDiameter(mesh)  # UFL symbolic representation of cell diameter
    h_func = dolfinx.fem.Function(Size)
    h_expr = dolfinx.fem.Expression(h, Size.element.interpolation_points())
    h_func.interpolate(h_expr)
    element_sizes = h_func.x.array
    filtered_element_ids = element_ids[element_sizes[element_ids] > mesh_size]
    return filtered_element_ids

L = ufl.dot(b,du)*ufl.dx + ufl.dot(t_bar,du)*ds_top
a = (degradation_fn(dold)+ k)*ufl.inner(sigmah(u),epsilon(du))*ufl.dx
L_frac = 2*dd*H(uold,unew,Hold) *ufl.dx
a_frac = G_cr*l*ufl.dot(ufl.grad(d),ufl.grad(dd)) *ufl.dx + ((G_cr/l) + 2*H(uold,unew,Hold))*ufl.inner(d,dd) * ufl.dx

problem1 = dolfinx.fem.petsc.LinearProblem(a, L, bcs = [bc_ubot, bc_utop, bc_utopy], petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type":"mumps"})
problem2 = dolfinx.fem.petsc.LinearProblem(a_frac, L_frac,bcs = [bc_dcrack], petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type":"mumps"})

load = 0#load = 0.3 #for initial u = 6e-3
incremental_tol=1e-4
trigger_tol=1e-3
max_iterations=50

beta1=1
en_threshold= (3/8)*(G_cr/l)
beta3 = 1.1  # Example scaling factor
radius = beta3 * l  # Maximum distance allowed
active_el_refine=2
added_el_refine=1
damaged_el_refine=2
alpha_threshold = 0.9
characteristic_element_size = l/3
#
while load <= 1:
    u_topval.value = load*total_u
    u_top.interpolate(u_topexpr)
    
    i = 1
    converged = False
    trigger = 0
    linesearch = False#True
    totenold=np.inf
    alpha_star=1
    flag_firstiter=0
    err = 1
        
    while i < max_iterations :
        
        flag_samemesh = False
        MPI.COMM_WORLD.Barrier()
        ###################################AMR LOOP######################################################
        if flag_firstiter == 0: 
            unew = problem1.solve()
            imesh =dmesh  #imesh,ict,ift = dmesh,ct,ft
            if load != 0:
                amesh = mesh
            V_en_elastic = fem.functionspace(dmesh, ("DG", 0))
            en_elastic_expr = fem.Expression(ufl.inner(sigmah(unew),epsilon(unew)), V_en_elastic.element.interpolation_points())
            en_elastic = fem.Function(V_en_elastic)
            en_elastic.interpolate(en_elastic_expr)
            energy_values = en_elastic.x.array
            active_el_cell_ids = np.flatnonzero(energy_values > en_threshold).astype(np.int32)
            
            if active_el_cell_ids.size > 0:
                midpoints_old=[]
                midpoints=[]
                for num in range(active_el_refine-added_el_refine):
                    #
                    if num==0:
                        mesh , element_list, midpoints = obtain_midpoints(dmesh) #imesh
                        active_el_cell_ids = filter_elements_by_size(mesh,active_el_cell_ids,characteristic_element_size)
                        mesh.topology.create_entities(1)
                        local_edges = dolfinx.mesh.compute_incident_entities(mesh.topology, active_el_cell_ids, mesh.topology.dim, 1)
                    else:
                        mesh , element_list, midpoints = obtain_midpoints(mesh)
                        mask = np.all(np.isin(midpoints, midpoints_old, assume_unique=False), axis=0)
                        refined_midpoints = midpoints[:, ~mask]
                        refined_cell_ids = element_list[~mask]
                        refined_cell_ids = filter_elements_by_size(mesh,refined_cell_ids,characteristic_element_size)
                        mesh.topology.create_entities(1)
                        local_edges = dolfinx.mesh.compute_incident_entities(mesh.topology, refined_cell_ids, mesh.topology.dim, 1)
                    mesh = dolfinx.mesh.refine(mesh, local_edges)
                    midpoints_old=midpoints
                #
                midpoints_old=[]
                midpoints=[]
                for num in range(added_el_refine):
                    if num==0:
                        mesh , element_list, midpoints = obtain_midpoints(mesh)
                        mesh_coarse , element_list_coarse, midpoints_coarse = obtain_midpoints(imesh)
                        tree = cKDTree(midpoints.T)  # midpoints is the full list of element midpoints
                        mask = np.all(np.isin(midpoints, midpoints_coarse, assume_unique=False), axis=0)
                        active_el_midpoints = midpoints[:, ~mask]
                        active_el_id = element_list[~mask]
                        neighbor_indices = tree.query_ball_point(active_el_midpoints.T, radius)
                        added_el_ids = np.unique(np.hstack(neighbor_indices))
                        added_el_ids = np.array(added_el_ids, dtype=np.int32)
                        added_el_ids = filter_elements_by_size(mesh,added_el_ids,characteristic_element_size)
                        mesh.topology.create_entities(1)
                        local_edges = dolfinx.mesh.compute_incident_entities(mesh.topology, added_el_ids, mesh.topology.dim, 1)
                    else:
                        mesh , element_list, midpoints = obtain_midpoints(mesh)
                        mask = np.all(np.isin(midpoints, midpoints_old, assume_unique=False), axis=0)
                        refined_midpoints = midpoints[:, ~mask]
                        refined_cell_ids = element_list[~mask]
                        refined_cell_ids = filter_elements_by_size(mesh,refined_cell_ids,characteristic_element_size)
                        mesh.topology.create_entities(1)
                        local_edges = dolfinx.mesh.compute_incident_entities(mesh.topology, refined_cell_ids, mesh.topology.dim, 1)
                    mesh = dolfinx.mesh.refine(mesh, local_edges)
                    midpoints_old=midpoints
            if active_el_cell_ids.size==0:
                mesh =dmesh
            ##alpha <alpha_s
            midpoints_old=[]
            midpoints=[]
            neighbor_indices=np.array([])
            for num in range(damaged_el_refine):
                #
                if num==0:
                    if load!=0:
                        amesh, aelement_list, amidpoints = obtain_midpoints(amesh)
                        d_phase = solr.sub(1).collapse()
                        d_values = d_phase.x.array
                        damaged_nodes = np.flatnonzero(d_values > alpha_threshold)
                        connectivity = amesh.topology.connectivity(0,amesh.topology.dim)
                        damaged_element_ids = set()
                        for node in damaged_nodes:
                            damaged_element_ids.update(connectivity.links(node))

                        damaged_element_ids = np.array(list(damaged_element_ids), dtype=np.int32)
                        damaged_midpoints=amidpoints[:, damaged_element_ids]                        
                        mesh , element_list, midpoints = obtain_midpoints(mesh)
                        tree = cKDTree(midpoints.T)  # KDTree for fast nearest neighbor search
                        imesh_size = np.min(element_size_list(mesh))  # Estimate element size
                        search_radius = imesh_size  # Adjust search radius based on refinement level
                        neighbor_indices = tree.query_ball_point(damaged_midpoints.T, search_radius)
                        if neighbor_indices.size == 0:
                            break
                        new_damaged_elements_ids = np.unique(np.hstack(neighbor_indices))
                        new_damaged_elements_ids=np.array(new_damaged_elements_ids,dtype=np.int32)
                        new_damaged_elements_ids = filter_elements_by_size(mesh,new_damaged_elements_ids,characteristic_element_size)
                        mesh.topology.create_entities(1)
                        local_edges = dolfinx.mesh.compute_incident_entities(mesh.topology, new_damaged_elements_ids, mesh.topology.dim, 1)
                    else:
                        mesh , element_list, midpoints = obtain_midpoints(mesh)
                        local_edges=dolfinx.mesh.locate_entities_boundary(mesh, 1, crack)
                else:
                    mesh , element_list, midpoints = obtain_midpoints(mesh)
                    mask = np.all(np.isin(midpoints, midpoints_old, assume_unique=False), axis=0)
                    refined_midpoints = midpoints[:, ~mask]
                    refined_cell_ids = element_list[~mask]
                    refined_cell_ids = filter_elements_by_size(mesh,refined_cell_ids,characteristic_element_size)
                    mesh.topology.create_entities(1)
                    local_edges = dolfinx.mesh.compute_incident_entities(mesh.topology, refined_cell_ids, mesh.topology.dim, 1)
                #
                mesh = dolfinx.mesh.refine(mesh, local_edges)
                midpoints_old=midpoints
        ##################################AMR LOOP####################################################
        if flag_firstiter == 0 :
            mix = basix.ufl.mixed_element([P,S])
            Vr = fem.functionspace(mesh,mix)
            Wr = fem.functionspace(mesh,R)
            dur , ddr = ufl.TestFunctions(Vr)
            solr = fem.Function(Vr)
            sololdr = fem.Function(Vr)
            HHr = fem.Function(Wr)
            Holdr = fem.Function(Wr) # history field at the previous time step
            ur, dr = ufl.split(solr)
            uoldr, doldr = ufl.split(sololdr)
            V0r,_ = Vr.sub(0).collapse()
            V01r, _ = Vr.sub(0).sub(0).collapse()
            V02r, _ = Vr.sub(0).sub(1).collapse()
            V1r,_ = Vr.sub(1).collapse()
            bot_facetsr = dolfinx.mesh.locate_entities_boundary(mesh, 1, bot)
            top_facetsr = dolfinx.mesh.locate_entities_boundary(mesh, 1, top)
            crack_facetsr = dolfinx.mesh.locate_entities_boundary(mesh, 1, crack)
            bot_dofsr = fem.locate_dofs_topological((Vr.sub(0),V0r), 1, bot_facetsr)
            top_dofsr = fem.locate_dofs_topological((Vr.sub(0).sub(0),V01r), 1, top_facetsr)
            top_dofs_yr = fem.locate_dofs_topological((Vr.sub(0).sub(1),V02r), 1, top_facetsr)
            crack_dofsr = fem.locate_dofs_topological((Vr.sub(1),V1r), 1, crack_facetsr)
            u_botr = fem.Function(V0r)
            u_botvalr = fem.Constant(mesh,ScalarType([0.0,0.0]))
            u_botexprr = fem.Expression(u_botvalr,V0r.element.interpolation_points())
            u_botr.interpolate(u_botexprr)
            u_topr = fem.Function(V01r)
            u_topvalr = fem.Constant(mesh,ScalarType(0.00))    
            u_topexprr = fem.Expression(u_topvalr,V01r.element.interpolation_points())
            u_topr.interpolate(u_topexprr)
            u_topyr = fem.Function(V02r)
            u_topyvalr = fem.Constant(mesh,ScalarType(0))    
            u_topyexprr = fem.Expression(u_topyvalr,V02r.element.interpolation_points())
            u_topyr.interpolate(u_topyexprr)
            d_crackr = fem.Function(V1r)
            d_crackr.x.array[:] = 1 
            bc_ubotr = fem.dirichletbc(u_botr,bot_dofsr,Vr.sub(0))
            bc_utopr = fem.dirichletbc(u_topr, top_dofsr, Vr.sub(0).sub(0))
            bc_utopyr = fem.dirichletbc(u_topyr, top_dofs_yr, Vr.sub(0).sub(1))
            bc_dcrackr = fem.dirichletbc(d_crackr, crack_dofsr, Vr.sub(1))
            t_barr = fem.Constant(mesh,np.array([0.0,0.0],dtype=ScalarType))
            br = fem.Constant(mesh,np.array([0.0,0.0],dtype=ScalarType))
            nr = ufl.FacetNormal(mesh)

            facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, top)
            values = np.full(len(facets), 2, dtype=np.int32)
            ft_manual = dolfinx.mesh.meshtags(mesh, 1, facets, values)
            ds_topr = ufl.Measure(integral_type='ds', domain=mesh, subdomain_data=ft_manual, subdomain_id=2) 

            mechr = ufl.inner(sigma(ur, doldr),epsilon(dur))*ufl.dx - ufl.dot(br,dur)*ufl.dx - ufl.dot(t_barr,dur)*ds_topr
            fracr = G_cr*l*ufl.dot(ufl.grad(dr),ufl.grad(ddr)) *ufl.dx + ((G_cr/l) + 2*H(uoldr,ur,Holdr))*ufl.inner(dr,ddr) * ufl.dx - 2*ddr*H(uoldr,ur,Holdr) *ufl.dx
            Fr = mechr + fracr

            sololdnrr = fem.Function(Vr) 
            uoldnrr, doldnrr = ufl.split(sololdnrr)
            solnrr = fem.Function(Vr) ###not used
            unrr, dnrr = ufl.split(solnrr)
            residualr = dolfinx.fem.form(Fr)
            dsolr = ufl.TrialFunction(Vr)
            Jr = ufl.derivative(Fr, solr, dsolr)
            jacobianr = dolfinx.fem.form(Jr)
            delsolr = fem.Function(Vr)
            delur, deldr = ufl.split(delsolr)
            delur = delsolr.sub(0)
            deldr = delsolr.sub(1)
            Ar = dolfinx.fem.petsc.create_matrix(jacobianr)
            Lr = dolfinx.fem.petsc.create_vector(residualr)
            solverr = PETSc.KSP().create(mesh.comm)
            solverr.setType(PETSc.KSP.Type.PREONLY)
            solverr.getPC().setType(PETSc.PC.Type.LU)
            optsr = PETSc.Options()
            prefixr = f"solver_{id(solverr)}"
            solverr.setOptionsPrefix(prefixr)
            option_prefixr = solverr.getOptionsPrefix()
            optsr[f"{option_prefixr}ksp_type"] = "preonly"
            optsr[f"{option_prefixr}pc_type"] = "lu"
            optsr[f"{option_prefixr}pc_factor_mat_solver_type"] = "mumps"
            solverr.setFromOptions()
            solverr.setOperators(Ar)
            ur = solr.sub(0)
            dr = solr.sub(1)
            uoldr = sololdr.sub(0)
            doldr = sololdr.sub(1)
            uoldnrr = sololdnrr.sub(0)
            doldnrr = sololdnrr.sub(1)
            unrr = solnrr.sub(0)
            dnrr = solnrr.sub(1)
        #
        MPI.COMM_WORLD.Barrier()
        u_topvalr.value = load*total_u
        u_topr.interpolate(u_topexprr)
        
        with Lr.localForm() as loc_L:
            loc_L.set(0)
        Ar.zeroEntries()
        dolfinx.fem.petsc.assemble_matrix(Ar, jacobianr, bcs=[bc_ubotr, bc_utopr, bc_dcrackr, bc_utopyr])
        Ar.assemble()
        dolfinx.fem.petsc.assemble_vector(Lr, residualr)
        Lr.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        Lr.scale(-1)
        dolfinx.fem.petsc.apply_lifting(Lr, [jacobianr], [[bc_ubotr, bc_utopr, bc_dcrackr, bc_utopyr]], x0=[solr.x.petsc_vec],scale=1)
        dolfinx.fem.petsc.set_bc(Lr, [bc_ubotr, bc_utopr, bc_dcrackr, bc_utopyr], solr.x.petsc_vec, 1.0)
        Lr.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
        solverr.solve(Lr, delsolr.x.petsc_vec)
        delsolr.x.scatter_forward()
        residual_normr = Lr.norm(1)
        correction_normr = delsolr.x.petsc_vec.norm(1)
        solr.x.array[:] += alpha_star*delsolr.x.array
        solr.x.scatter_forward()
        sololdnrr.x.array[:] = solr.x.array  ##not used. can remove
        sololdnrr.x.scatter_forward()
        i += 1
        flag_firstiter=1      
        
        if residual_normr < incremental_tol or correction_normr < incremental_tol: 
            converged=True
            break

    H_exprr = fem.Expression(ufl.conditional(ufl.operators.gt(psi_e(ur),Holdr),psi_e(ur),Holdr),Wr.element.interpolation_points())
    Holdr.interpolate(H_exprr)
    Holdr.x.scatter_forward()
    
    sololdr.x.array[:] = solr.x.array  
    sololdr.x.scatter_forward()
    solold_carry = sololdr.copy()  
    
    filename = f"output_step_{load*total_u}.xdmf"
    with io.XDMFFile(MPI.COMM_WORLD, filename, "w") as xdmf:
        mesh.name = "mesh"
        xdmf.write_mesh(mesh)
        xdmf.write_function(solr.sub(1).collapse(), t=load*total_u)
    
    if not converged:
        if MPI.COMM_WORLD.rank == 0:
            print("Solver did not converge! LOL.",flush=True)
        break
    load += delta_u
    #
# end stopwatch
MPI.COMM_WORLD.Barrier()
end_time = MPI.Wtime()
print(f"Elapsed time: {end_time-start_time} s.",flush=True) #remove later

I think all problems are within AMR algorithm and improper communication between processes for mesh.

And here is the gmsh file

SetFactory("OpenCASCADE");
Point(1) = {0, 0, 0, 1};
Point(2) = {0, 1, 0, 1};
Point(3) = {1, 1, 0, 1};
Point(4) = {1, 0, 0, 1};
Point(5) = {0.5, 0.5, 0, 1};
Point(6) = {0, 0.5+ 0.001, 0, 1};
Point(7) = {0, 0.5-0.001, 0, 1};

Line(1) = {1, 4};
Line(2) = {4, 3};
Line(3) = {3, 2};
Line(4) = {2, 6};
Line(5) = {6, 5};
Line(6) = {5, 7};
Line(7) = {7, 1};
Curve Loop(1) = {1,2,3,4,5,6,7}
Plane Surface(1) = {1};

Physical Curve("top", 2) = {3};
Physical Curve("bottom", 1) = {1};
Physical Curve("crack_top", 31) = {5};
Physical Curve("crack_bot", 32) = {6};
Physical Surface("surf", 4) = {1};

Field[1] = Box;
Field[1].Thickness = 0.2;
Field[1].VIn = 0.025;
Field[1].VOut = 0.025;
Field[1].XMax = 1;
Field[1].XMin = 0;
Field[1].YMax = 1;
Field[1].YMin = 0;
Background Field = 1;
Field[3] = Min;
Field[3].FieldsList = {1};
Background Field = 3;
//2d algo- frontdelau for quad (eexperimental)
//recomined algo-simple fullquad
//rrecombine all triangular
//+

In serial, after the first load step, only the pre-crack edges are refined, which is correct:

However, in parallel(here 3 cores), the refinement behaves unexpectedly, even at the first step.

can someone suggest how I should go about diagnosing and resolving this? Any guidance would be greatly appreciated!

Thanks in advance for your help.

This example is unfortunately very long and difficult for volunteers on this forum to parse quickly. If it could be condensed into the salient components of the issue in a minimal working example, you’d be more likely to receive advice.