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.