# Issue with Mesh Study in Elasticity Problem: Inconsistent Results with Increased Elements

Hey all, I’m working on a simple elasticity problem and encountering an issue during my mesh study. The code runs fine with a low number of elements (39904), but when I increase the number of elements (158564 elements), I get constant results across the entire domain. Here are the two meshes I’m using

The first mesh works well but the second mesh doesn’t.

Here are the outputs for both meshes:

I don’t know what the problem is. Is that because the direct solver that I am using?

And here is the code:

``````from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista
import math
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological, assemble, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary
from petsc4py.PETSc import ScalarType
import meshio
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, split)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type

# .........................................  Non-Dimensionalization Parameters
Fl = 1e-3
Fs = 1e-5
Fp = 1e-4
# # # ....................
C_star = Fl ** 2 / Fs
h_star = 1 / Fl
p_star = Fl ** 3
F_star = Fs
# -------------------------Material Properties
E = 1 * 130e9 * C_star
nu = 0.28
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
h = 0.76e-3 / Fl
d1 = 1.13e-3 / Fl
tb = 0.1e-3 / Fl
########################
H = h + tb
mesh, cell_markers, facet_markers = gmshio.read_from_msh("Pyramid.msh", MPI.COMM_WORLD,
gdim=3)  # gmshio.read_from_msh always needs physical group in .msh file
Force = 200.00000 / Fs
Area2 = 1.2769000000000007e-06 / Fl ** 2
tr_S2 = Constant(mesh, ScalarType((0, 0, -Force / Area2)))
n = FacetNormal(mesh)
MIXEL = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))  # Displacement Element
Space = functionspace(mesh, MIXEL)
# ......................................... Define Trail and Test Functions
u = TrialFunction(Space)
del_u = TestFunction(Space)
# ....................................... Mark boundary regions
tol = 1 * 1e-16

def Bottom(x):
return np.isclose(x[2], 0, atol=tol)

# ....................................... ds
boundaries = [(1, lambda x: np.isclose(x[2], 0.0, atol=tol)), (2, lambda x: np.isclose(x[2], H, atol=tol))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

# # ...................................... BCs
# Clamp BC for left side
clamp_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Bottom)
clamp_dofs = locate_dofs_topological(Space, mesh.topology.dim - 1, clamp_facets)
bc0 = dirichletbc(np.array([0, 0, 0], dtype=default_scalar_type), clamp_dofs, Space)
bc = [bc0]  # BC 4
# ......................................... Define tensor
delta = Identity(3)
i, j, k, l = indices(4)

# Strain Tensor
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])

# Stress Tensor
def Elas(la, mu):
return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l],
(i, j, k, l))  # Elasticity tensor

def StressT(u):
return as_tensor((Elas(la, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])

# ......................................... Define Variational Formulation
VarForm = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
a = VarForm
L = tr_S2[i] * del_u[i] * ds(2)
bilinear_form = form(a)
linear_form = form(L)
A = assemble_matrix(bilinear_form, bcs=bc)
A.assemble()
b = assemble_vector(linear_form)
apply_lifting(b, [bilinear_form], [bc])
set_bc(b, bc)

problem = LinearProblem(a, L, bcs=bc,
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
Result = problem.solve()

with VTXWriter(mesh.comm, ("Pyramid/DIS_Pyramid.bp"), [Result]) as vtx:
vtx.write(0.0)

``````

That’s possible. Please try using `superlu` or `superlu_dist` and check if you still get the same result.

Thanks for the reply, I forgot to provide the mesh. The following code generates the mesh in GMSH.

``````
import gmsh
import sys

##################
gmsh.initialize()
# =====================================    Goccmetry      ======================
import gmsh
import sys
import math

gmsh.initialize()
Fl = 1e-3
lc = 0.0
h = 0.76e-3 / Fl
d1 = 1.13e-3 / Fl
d2 = d1
lb = 0.1e-3 / Fl
tb = 0.1e-3 / Fl
theta = math.atan(h / (0.5 * (0.00272 / Fl - d1)))  #54.74
Points = [[0, 0, 0], [2 * lb + d1 + 2 * h / math.tan(theta), 0, 0],
[2 * lb + d1 + 2 * h / math.tan(theta), 2 * lb + d2 + 2 * h / math.tan(theta), 0],
[0, 2 * lb + d2 + 2 * h / math.tan(theta), 0]]
print(d1 + 2 * h / math.tan(theta))
PN = 1  # point Number
LN = 1  # Line Number
SN = 1  # Surface Number
VN = 1  # Volum Number
Vol_list = []
Sur_list = []
Sur_list2 = []
Sur_list3 = []
Sur_list4 = []
Sur_list5 = []
Sur_list6 = []
Edge_list_1 = []  # top edges
Edge_list_2 = []  # Bottom edges
Edge_list_3 = []  # Side edges
Edge_list_4 = []  # in-plane edges
Edge_list_5 = []  # corner edges
Edge_list_6 = []  # Base edges
i=0
gmsh.model.geo.addPoint(Points[i][0] + lb + h / math.tan(theta), Points[i][1] + lb + h / math.tan(theta), tb + h,
lc, PN)
PN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + h / math.tan(theta) + d1, Points[i][1] + lb + h / math.tan(theta),
tb + h, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
Edge_list_1.append(LN)
LN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + h / math.tan(theta) + d1, Points[i][1] + lb + h / math.tan(theta) + d2,
tb + h, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
Edge_list_1.append(LN)
LN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + h / math.tan(theta), Points[i][1] + lb + h / math.tan(theta) + d2,
tb + h, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
Edge_list_1.append(LN)
LN += 1
gmsh.model.geo.addLine(PN - 1, PN - 4, LN)
Edge_list_1.append(LN)
LN += 1
#########
gmsh.model.geo.addPoint(Points[i][0] + lb, Points[i][1] + lb, tb, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + 2 * h / math.tan(theta) + d1, Points[i][1] + lb, tb, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
Edge_list_2.append(LN)
LN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + 2 * h / math.tan(theta) + d1,
Points[i][1] + lb + 2 * h / math.tan(theta) + d2, tb, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
Edge_list_2.append(LN)
LN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb, Points[i][1] + lb + 2 * h / math.tan(theta) + d2, tb, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
Edge_list_2.append(LN)
LN += 1
gmsh.model.geo.addLine(PN - 1, PN - 4, LN)
Edge_list_2.append(LN)
LN += 1
####################
gmsh.model.geo.addLine(PN - 8, PN - 4, LN)
Edge_list_3.append(LN)
LN += 1
gmsh.model.geo.addLine(PN - 7, PN - 3, LN)
Edge_list_3.append(LN)
LN += 1
gmsh.model.geo.addLine(PN - 6, PN - 2, LN)
Edge_list_3.append(LN)
LN += 1
gmsh.model.geo.addLine(PN - 5, PN - 1, LN)
Edge_list_3.append(LN)
LN += 1
####################
gmsh.model.geo.addCurveLoop([LN - 12, LN - 11, LN - 10, LN - 9], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 12, (LN - 3), -(LN - 8), -(LN - 4)], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 11, (LN - 2), -(LN - 7), -(LN - 3)], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 10, (LN - 1), -(LN - 6), -(LN - 2)], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 9, (LN - 4), -(LN - 5), -(LN - 1)], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 8, LN - 7, LN - 6, LN - 5], SN)
SN += 1
Sur_list.append(SN - 1)
print(SN)
gmsh.model.geo.addSurfaceLoop([SN - 1, SN - 2, SN - 3, SN - 4, SN - 5, SN - 6], VN)
Vol_list.append(VN)
VN += 1
PN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + 2 * h / math.tan(theta) + d1 + lb, Points[i][1], tb, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + 2 * h / math.tan(theta) + d1 + lb,
Points[i][1] + 2 * lb + 2 * h / math.tan(theta) + d2, tb, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Points[i][0], Points[i][1] + 2 * lb + 2 * h / math.tan(theta) + d2, tb, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 8, PN - 4, LN)
LN += 1
Edge_list_4.append(LN - 1)
gmsh.model.geo.addLine(PN - 7, PN - 3, LN)
LN += 1
Edge_list_4.append(LN - 1)
gmsh.model.geo.addLine(PN - 6, PN - 2, LN)
LN += 1
Edge_list_4.append(LN - 1)
gmsh.model.geo.addLine(PN - 5, PN - 1, LN)
LN += 1
Edge_list_4.append(LN - 1)
gmsh.model.geo.addLine(PN - 4, PN - 3, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 3, PN - 2, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 1, PN - 4, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addCurveLoop([LN - 16, LN - 7, -(LN - 4), -(LN - 8)], SN)
SN += 1
Sur_list2.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 15, LN - 6, -(LN - 3), -(LN - 7)], SN)
SN += 1
Sur_list3.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 14, LN - 5, -(LN - 2), -(LN - 6)], SN)
SN += 1
Sur_list4.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 13, LN - 8, -(LN - 1), -(LN - 5)], SN)
SN += 1
Sur_list5.append(SN - 1)

gmsh.model.geo.synchronize()
Corner_points = [PN - 4, PN - 3, PN - 2, PN - 1]
Left_down_corner = gmsh.model.getValue(0, Corner_points[0], [])
Right_down_corner = gmsh.model.getValue(0, Corner_points[1], [])
Right_up_corner = gmsh.model.getValue(0, Corner_points[2], [])
Left_up_corner = gmsh.model.getValue(0, Corner_points[3], [])

PN += 1
PN += 1
PN += 1
PN += 1
gmsh.model.geo.addLine(PN - 8, PN - 4, LN)
LN += 1
Edge_list_6.append(LN-1)
gmsh.model.geo.addLine(PN - 7, PN - 3, LN)
LN += 1
Edge_list_6.append(LN - 1)
gmsh.model.geo.addLine(PN - 6, PN - 2, LN)
LN += 1
Edge_list_6.append(LN - 1)
gmsh.model.geo.addLine(PN - 5, PN - 1, LN)
LN += 1
Edge_list_6.append(LN - 1)
gmsh.model.geo.addLine(PN - 4, PN - 3, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 3, PN - 2, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 1, PN - 4, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addCurveLoop([LN - 12, LN - 7, -(LN - 4), -(LN - 8)], SN)
SN += 1
Sur_list2.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 11, LN - 6, -(LN - 3), -(LN - 7)], SN)
SN += 1
Sur_list3.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 10, LN - 5, -(LN - 2), -(LN - 6)], SN)
SN += 1
Sur_list4.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 9, LN - 8, -(LN - 1), -(LN - 5)], SN)
SN += 1
Sur_list5.append(SN - 1)
gmsh.model.geo.synchronize()
Corner_points = [PN - 12, PN - 11, PN - 10, PN - 9]
Left_down_corner = gmsh.model.getValue(0, Corner_points[0], [])
Right_down_corner = gmsh.model.getValue(0, Corner_points[1], [])
Right_up_corner = gmsh.model.getValue(0, Corner_points[2], [])
Left_up_corner = gmsh.model.getValue(0, Corner_points[3], [])
PN += 1
PN += 1
PN += 1
PN += 1
gmsh.model.geo.addLine(PN - 8, PN - 4, LN)
LN += 1
Edge_list_4.append(LN-1)
gmsh.model.geo.addLine(PN - 7, PN - 3, LN)
LN += 1
Edge_list_4.append(LN - 1)
gmsh.model.geo.addLine(PN - 6, PN - 2, LN)
LN += 1
Edge_list_4.append(LN - 1)
gmsh.model.geo.addLine(PN - 5, PN - 1, LN)
LN += 1
Edge_list_4.append(LN - 1)
# ################3
gmsh.model.geo.addLine(PN - 4, PN - 16, LN)
LN += 1
Edge_list_6.append(LN - 1)
gmsh.model.geo.addLine(PN - 3, PN - 15, LN)
LN += 1
Edge_list_6.append(LN - 1)
gmsh.model.geo.addLine(PN - 2, PN - 14, LN)
LN += 1
Edge_list_6.append(LN - 1)
gmsh.model.geo.addLine(PN - 1, PN - 13, LN)
LN += 1
Edge_list_6.append(LN - 1)
################
gmsh.model.geo.addLine(PN - 4, PN - 3, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 3, PN - 2, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 1, PN - 4, LN)
LN += 1
Edge_list_5.append(LN - 1)
# ################################
gmsh.model.geo.addCurveLoop([LN - 28, LN - 20, (LN - 12), (LN - 8)], SN)
SN += 1
Sur_list2.append(SN - 1)
Sur_list5.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 16, LN - 11, -(LN - 4), -(LN - 12)], SN)
SN += 1
Sur_list2.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 27, LN - 19, (LN - 11), (LN - 7)], SN)
SN += 1
Sur_list2.append(SN - 1)
Sur_list3.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 4, LN - 7, -(LN - 36), -(LN - 8)], SN)
SN += 1
Sur_list2.append(SN - 1)

Vol_list.append(VN)
VN += 1
############################
gmsh.model.geo.addCurveLoop([LN - 3, LN - 6, -(LN - 35), -(LN - 7)], SN)
SN += 1
Sur_list3.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 15, LN - 10, -(LN - 3), -(LN - 11)], SN)
SN += 1
Sur_list3.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 10, LN - 6, (LN - 26), (LN - 18)], SN)
SN += 1
Sur_list3.append(SN - 1)
Vol_list.append(VN)
VN += 1
Sur_list4.append(SN - 1)
################################
gmsh.model.geo.addCurveLoop([LN - 14, LN - 9, -(LN - 2), -(LN - 10)], SN)
SN += 1
Sur_list4.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 2, LN - 5, -(LN - 34), -(LN - 6)], SN)
SN += 1
Sur_list4.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 9, LN - 5, (LN - 25), (LN - 17)], SN)
SN += 1
Sur_list4.append(SN - 1)
Vol_list.append(VN)
VN += 1
Sur_list5.append(SN - 1)
###############################
gmsh.model.geo.addCurveLoop([LN - 13, LN - 12, -(LN - 1), -(LN - 9)], SN)
SN += 1
Sur_list5.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 1, LN - 8, -(LN - 33), -(LN - 5)], SN)
SN += 1
Sur_list5.append(SN - 1)

Vol_list.append(VN)
VN += 1
######################################
Sur_list6 = [SN-21, SN-9, SN-8, SN-4, SN-1]
gmsh.model.geo.addCurveLoop([LN - 4, LN - 3, (LN - 2), (LN - 1)], SN)
SN += 1
Sur_list6.append(SN - 1)

Vol_list.append(VN)
VN += 1
# =====================================    Mesh Option      ======================
## UP
N_up = 50
type_up = "Bump"
R_up = 0.5
## Bottom
N_bottom = N_up
type_bottom = "Bump"
R_bottom = 0.5
## Side
N_side = 60
type_side = "Bump"
R_side = 1
####
N_inPlane = 6
###
N_base= 6
# ----------------------- UP
for i in Edge_list_1:
gmsh.model.geo.mesh.setTransfiniteCurve(i, N_up, type_up, R_up)
for i in Edge_list_2:
gmsh.model.geo.mesh.setTransfiniteCurve(i, N_bottom, type_bottom, R_bottom)
for i in Edge_list_3:
gmsh.model.geo.mesh.setTransfiniteCurve(i, N_side, type_side, R_side)

for i in Edge_list_4:
gmsh.model.geo.mesh.setTransfiniteCurve(i, N_inPlane)#, type_side, R_side)

for i in Edge_list_5:
gmsh.model.geo.mesh.setTransfiniteCurve(i, N_bottom)#, type_side, R_side)

for i in Edge_list_6:
gmsh.model.geo.mesh.setTransfiniteCurve(i, N_base)#, type_side, R_side)

for i in Sur_list + Sur_list2 +Sur_list3 + Sur_list4 + Sur_list5 + Sur_list6:
gmsh.model.geo.mesh.setTransfiniteSurface(i)

for i in Sur_list + Sur_list2 +Sur_list3 + Sur_list4 + Sur_list5 + Sur_list6:
gmsh.model.geo.mesh.setRecombine(2, i)

for i in Vol_list:
gmsh.model.geo.mesh.setTransfiniteVolume(i)

gmsh.model.geo.remove_all_duplicates()
gmsh.model.geo.synchronize()
gmsh.option.setNumber('Mesh.MeshSizeMax', 0.005)
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(3)
# ##############
gmsh.write("Pyramid.msh")
######################
if '-nopopup' not in sys.argv:
gmsh.fltk.run()
#################################
gmsh.finalize()

``````

@francesco-ballarin , When I use the `superlu_dist` the code terminates due to the low memory error, even with the first mesh that has lower number of elements. And when I use `superlu`, I get the following

Thanks

Then I would guess you are running out of RAM. Is there any way you can try the same code on a workstation with more RAM?

Thanks @francesco-ballarin for your reply. I can’t upgrade the RAM on my PC at the moment, but I can plan to get a system with higher RAM in the future. However, the issue is that it doesn’t work for a mesh with a lower number of elements (39904), which I can run with `mumps`. Does this mean that using `superlu_dist` or `superlu` works better for large systems (high DOF) but requires more RAM than using `mumps`?

Do you think there’s no problem with my code, and that the strange results in my sumulation are just something that happens in large models when using `mumps`?

I have one more question. I’m currently using a direct solver in this code. Should I try using an iterative solver instead? Could that potentially solve the issue?

Thanks,

If the issue is the direct solver and/or the RAM, then using an iterative solver could help. Unfortunately I can’t help much further, since neither solid mechanics nor iterative solvers are my field of expertise. You may want to wait for further comments by someone more knowledgeable than me in those fields.

Thanks @francesco-ballarin for your time. I will wait to see if I get more comments.

Hey all,
I still need help with this issue and I am not sure if going after setting up an iterative solver is the best way to solve this problem. I would appreciate any comments on that.
Thanks,

Could you print the output of `problem.solver.getConvergedReason()` after the linear solver has given you a constant solution.

As a follow-up, how much RAM and how many processors do you have, and how many processes are you using with mpirun?

the out put is : 4
This is a linear problem. Shouldn’t that be 1?

I have 128 GB RAM and If i use one processor it is much faster (I don’t know why), but I can use up to 8.

Why is that super slow when I use more processors?

And here is the out put I get from
`dolfinx.common.list_timings(mesh.comm,[dolfinx.common.TimingType.wall]) `

Thanks

I do not have the sufficient memory on my laptop to run the code, so this would have to wait for me until I am on a bigger system.

How did you install dolfinx? using conda, docker, spack or something else?

Btw. your problem has `3 922 347` dofs, which is really quite alot to throw at a direct solver.

I installed it using conda. And, currently I am working on setting an iterative solver for my problem.

Dear Dokken,

I decided to implement an iterative solver instead of a direct solver.

I followed one of your tutorials, here , and successfully ran my code. However, I encountered a couple of issues :

1- With the same mesh, the code is significantly slower with the iterative solver, taking about 25,580 seconds compared to 492 seconds with the direct solver.

2- The results from the iterative solver are worse than those from the direct solver. Using the same mesh, the outcomes are not the same.

Could you please provide some guidance on how to address these issues? Thank you for your assistance.

Here is the code I am running:

``````#  importing modules
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista
import math
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological, assemble, locate_dofs_topological, extract_function_spaces,
bcs_by_block)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem, \
assemble_matrix_nest, assemble_vector_nest, apply_lifting_nest, set_bc_nest
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_box, locate_entities, meshtags, locate_entities_boundary, CellType
# from petsc4py.PETSc import ScalarType
import meshio
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, split, action,
SpatialCoordinate)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type
import dolfinx

# ......... Define parameteres
Fl = 1e-8
Fs = 1e-5
Fp = 1e-4
# ....................
C_star = Fl ** 2 / Fs
h_star = 1 / Fl
mu_star = Fp / Fs
perm_star = Fp ** 2 / Fs
F_star = Fs
# -------------------------Material Properties
# ...............Silicon
E = 1 * 152e9 * C_star
nu = 0.33
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
perm = 146.6e-9 * perm_star
epsilon = 11e-9
non_local = 10e-9 / Fl
mu11 = 0.5 * 121e-6 * mu_star
mu12 = 0.5 * 121e-6 * mu_star
mu44 = 0.0 * mu_star
h = 0.76e-3 / Fl
H = h
mesh = create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([h, h, h])],
[20, 20, 20], cell_type=CellType.hexahedron)
# ....................................... Mark boundary regions
tol = 1e-16

def Bottom(x):
return np.isclose(x[2], 0, atol=tol)

def Top(x):
return np.isclose(x[2], H, atol=tol)

boundaries = [(1, lambda x: np.isclose(x[2], 0.0, atol=tol)), (2, lambda x: np.isclose(x[2], H, atol=tol))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
Force = 200.00000 / Fs
Area1 = 7.398400000000025e-06 / Fl ** 2  # Bottom
Area2 = 1.2769000000000007e-06 / Fl ** 2
tr_S1 = Constant(mesh, default_scalar_type((0, 0, Force / Area1)))  # Bottom
tr_S2 = Constant(mesh, default_scalar_type((0, 0, -Force / Area2)))
noder = FacetNormal(mesh)
# ------------------------------  Defining the variational form
# .... Define Elements
u_el = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))  # Displacement Element
p_el = element("Lagrange", mesh.basix_cell(), 1)  # Electric Potential Element
s_el = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3))  # Relaxed strain Element
lm_el = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3))  # Lagrange multipliers Element
# .... Define Function Spaces
V = functionspace(mesh, u_el)
Q = functionspace(mesh, p_el)
VG = functionspace(mesh, s_el)
AL = functionspace(mesh, lm_el)
# .... Define Trail and Test Functions
u, p, gradu, lamu = TrialFunction(V), TrialFunction(Q), TrialFunction(VG), TrialFunction(AL)
del_u, del_p, del_gradu, del_lamu = TestFunction(V), TestFunction(Q), TestFunction(VG), TestFunction(AL)
# ..... Create  forms
delta = Identity(len(u))
i, j, k, l, m, n = indices(6)
fl1 = np.zeros((3, 3, 3, 3))
fl1[0, 0, 0, 0] = 1
fl1[1, 1, 1, 1] = 1
fl1[2, 2, 2, 2] = 1
flexo1 = Constant(mesh, PETSc.ScalarType(fl1))
flexo = as_tensor(((mu11 - mu12 - 2 * mu44) * flexo1[i, j, k, l] + mu12 * delta[i, j] * delta[
k, l] + mu12 * delta[i, k] * delta[j, l] + mu12 * delta[i, l] * delta[j, k]), (i, j, k, l))

# .....   bilinear_form
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])

# Stress Tensor
def Elas(la, mu):
return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l],
(i, j, k, l))  # Elasticity tensor

def StressT(u):
return as_tensor((Elas(la, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])

return as_tensor((1. / 2. * (gradu[k, j].dx(i) + gradu[j, k].dx(i))), [i, j, k])

def SGElas(la, non_local, mu):
return as_tensor((la * (non_local ** 2) * delta[i, j] * delta[l, m] * delta[k, n] +
2 * mu * (non_local ** 2) * delta[i, l] * delta[j, m] * delta[k, n]),
(i, j, k, l, m, n))

kuu = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
kua = -lamu[j, k] * del_u[k].dx(j) * dx
kpp = -perm * delta[i, j] * (-1 * nabla_grad(p))[j] * (-1 * nabla_grad(del_p)[i]) * dx
ksp = -1 * flexo[l, i, j, k] * (-1 * nabla_grad(p))[l] * RStrianGrad(del_gradu)[i, j, k] * dx
ksa = lamu[j, k] * del_gradu[j, k] * dx
kau = del_lamu[j, k] * u[k].dx(j) * dx
kas = del_lamu[j, k] * gradu[j, k] * dx
Zer_AL = Constant(AL.mesh, default_scalar_type(0.0))
kaa = Zer_AL * lamu[j, k] * del_lamu[j, k] * dx
a = form([[kuu, None, None, kua], [None, kpp, kps, None], [None, ksp, kss, ksa], [kau, None, kas, kaa]])
# ..... create linear_form
ZerVG = Constant(VG.mesh, default_scalar_type((0, 0, 0)))
Zer_AL2 = Constant(AL.mesh, default_scalar_type((0, 0, 0)))
fu = tr_S2[i] * del_u[i] * ds(2)
fp = inner(Constant(Q.mesh, 0.), del_p) * dx
fs = ZerVG[i] * del_gradu[i, j] * noder[j] * ds
fa = Zer_AL2[i] * del_lamu[i, j] * noder[j] * ds
L = form([fu, fp, fs, fa])
# ....................................... Create BCs
Bottom_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Bottom)
# ............................
clamp = Constant(V.mesh, [0., 0., 0.])
clamp_dofs = locate_dofs_topological(V, mesh.topology.dim - 1, Bottom_facets)
bc0 = dirichletbc(clamp, clamp_dofs, V)
# .........Ground BC
ground = Constant(Q.mesh, 0.0)
ground_dofs = locate_dofs_topological(Q, mesh.topology.dim - 1, Bottom_facets)
bc1 = dirichletbc(ground, ground_dofs, Q)
# ....
bcs = [bc0, bc1]
# ------------------------------  Assemble the nested system
A = assemble_matrix_nest(a, bcs=bcs)
A.assemble()
b = assemble_vector_nest(L)
apply_lifting_nest(b, a, bcs=bcs)
spaces = extract_function_spaces(L)
bcs0 = bcs_by_block(spaces, bcs)
set_bc_nest(b, bcs0)
# ------------------------------  Create Subspace solver
ksp = PETSc.KSP().create(MPI.COMM_WORLD)
ksp.setOperators(A)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
###  SOLVING
Displ, potent, DisGrad, LaMA = Function(V), Function(Q), Function(VG), Function(AL)
Result = PETSc.Vec().createNest([Displ.vector, potent.vector, DisGrad.vector, LaMA.vector ])
ksp.solve(b, Result)
#...........
Displ.x.scatter_forward()
potent.x.scatter_forward()
LaMA.x.scatter_forward()

with VTXWriter(mesh.comm, ("DIS.bp"), [Displ]) as vtx:
vtx.write(0.0)

with VTXWriter(mesh.comm, ("POTbp"), [potent]) as vtx:
vtx.write(0.0)

``````

Here is one of the out put for

direct

and iterative solver

As I cannot run your problem on my default systems due to memory limitations, it’s hard to give any further guidance.

Have you check that either solver actually converges?
I.e. what is the output of: `ksp.getConvergedReason()`

Thanks for the reply. I think it is possible to run the code on your system. I run it with the following mesh

``````mesh = create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([h, h, h])],
[6, 6, 6], cell_type=CellType.hexahedron)
``````

and it takes only 26 seconds to complete. And here is the output of `ksp.getConvergedReason()`:

Direct Solver: `4`
Iterative Solver: `-3`

That means that the iterative solver doesn’t converge within the set number of iterations:
https://petsc.org/release/manualpages/KSP/KSPConvergedReason/

Since you are using nest matrices, you should really use a fieldsplit if you want an iterative solver:
https://petsc.org/release/manualpages/PC/PCFIELDSPLIT/
See for instance:
https://docs.fenicsproject.org/dolfinx/v0.8.0/python/demos/demo_stokes.html#nested-matrix-solver
for an example.

Thanks for the reply. It was helpful and I was able to edit my code based on the material you suggested. Now the output for `ksp.getConvergedReason()` is 2.

The code solve more a less correctly for the first unknow (Displacment) but it generated random and meaningless results for the other unknow(Potential). It seems to me that the iterative solver works fine only for the first block. Here is the output for the potential.

Do you have any idea what might be the problem and what I should work on?

Here is the updated code:

``````#  importing modules
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista
import math
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological, assemble, locate_dofs_topological, extract_function_spaces,
bcs_by_block)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector_nest, set_bc, LinearProblem, \
assemble_matrix_nest, assemble_vector_nest, apply_lifting_nest, set_bc_nest
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_box, locate_entities, meshtags, locate_entities_boundary, CellType
# from petsc4py.PETSc import ScalarType
import meshio
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, split, action,
SpatialCoordinate)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type
import time
import dolfinx
# ......... Define parameteres
Fl = 1e-8
Fs = 1e-5
Fp = 1e-4
# ....................
C_star = Fl ** 2 / Fs
h_star = 1 / Fl
mu_star = Fp / Fs
perm_star = Fp ** 2 / Fs
q_star = Fp / (Fs * Fl)
mo_star = Fp / (Fl ** 2)
diff_star = 1 / (Fl ** 2)
p_star = Fl ** 3
F_star = Fs
qe = (1.6021 * 1e-19) * q_star
# -------------------------Material Properties
# ...............Silicon
E = 1 * 152e9 * C_star
nu = 0.33
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
perm = 146.6e-9 * perm_star
epsilon = 11e-9
non_local = 10e-20 / Fl
mu11 = 0.5 * 121e-6 * mu_star
mu12 = 0.5 * 121e-6 * mu_star
mu44 = 0.0 * mu_star
h = 0.76e-3 / Fl
H = h
mesh = create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([h, h, h])],
[10, 10, 10], cell_type=CellType.hexahedron)
# ....................................... Mark boundary regions
tol = 1e-16
def Bottom(x):
return np.isclose(x[2], 0, atol=tol)

def Top(x):
return np.isclose(x[2], H, atol=tol)

boundaries = [(1, lambda x: np.isclose(x[2], 0.0, atol=tol)), (2, lambda x: np.isclose(x[2], H, atol=tol))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
Force = 200.00000 / Fs
Area1 = 7.398400000000025e-06 / Fl ** 2  # Bottom
Area2 = 1.2769000000000007e-06 / Fl ** 2
tr_S1 = Constant(mesh, default_scalar_type((0, 0, Force / Area1)))  # Bottom
tr_S2 = Constant(mesh, default_scalar_type((0, 0, -Force / Area2)))
noder = FacetNormal(mesh)
# ------------------------------  Defining the variational form
# .... Define Elements
u_el = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))  # Displacement Element
p_el = element("Lagrange", mesh.basix_cell(), 1)  # Electric Potential Element
s_el = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3))  # Relaxed strain Element
lm_el = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3))  # Lagrange multipliers Element
# .... Define Function Spaces
V = functionspace(mesh, u_el)
Q = functionspace(mesh, p_el)
VG = functionspace(mesh, s_el)
AL = functionspace(mesh, lm_el)
# .... Define Trail and Test Functions
u, p, gradu, lamu = TrialFunction(V), TrialFunction(Q), TrialFunction(VG), TrialFunction(AL)
del_u, del_p, del_gradu, del_lamu = TestFunction(V), TestFunction(Q), TestFunction(VG), TestFunction(AL)
# ..... Create  forms
delta = Identity(len(u))
i, j, k, l, m, n = indices(6)
# Flexoelectricity Tensor
fl1 = np.zeros((3, 3, 3, 3))
fl1[0, 0, 0, 0] = 1
fl1[1, 1, 1, 1] = 1
fl1[2, 2, 2, 2] = 1
flexo1 = Constant(mesh, PETSc.ScalarType(fl1))
flexo = as_tensor(((mu11 - mu12 - 2 * mu44) * flexo1[i, j, k, l] + mu12 * delta[i, j] * delta[
k, l] + mu12 * delta[i, k] * delta[j, l] + mu12 * delta[i, l] * delta[j, k]), (i, j, k, l))

# .....   bilinear_form
# Strain Tensor
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])

# Stress Tensor
def Elas(la, mu):
return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l],
(i, j, k, l))  # Elasticity tensor

def StressT(u):
return as_tensor((Elas(la, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])

# Strain Gradient Tensor using relaxed strain
return as_tensor((1. / 2. * (gradu[k, j].dx(i) + gradu[j, k].dx(i))), [i, j, k])

# Higher Order Stress
def SGElas(la, non_local, mu):
return as_tensor((la * (non_local ** 2) * delta[i, j] * delta[l, m] * delta[k, n] +
2 * mu * (non_local ** 2) * delta[i, l] * delta[j, m] * delta[k, n]),
(i, j, k, l, m, n))  # Strain Gradient Elasticity Tensor

kuu = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
kua = -lamu[j, k] * del_u[k].dx(j) * dx
kpp = -perm * delta[i, j] * (-1 * nabla_grad(p))[j] * (-1 * nabla_grad(del_p)[i]) * dx
ksp = -1 * flexo[l, i, j, k] * (-1 * nabla_grad(p))[l] * RStrianGrad(del_gradu)[i, j, k] * dx
ksa = lamu[j, k] * del_gradu[j, k] * dx
kau = del_lamu[j, k] * u[k].dx(j) * dx
kas = del_lamu[j, k] * gradu[j, k] * dx
Zer_AL = Constant(AL.mesh, default_scalar_type(0.0))
kaa = Zer_AL * lamu[j, k] * del_lamu[j, k] * dx
a = form([[kuu, None, None, kua], [None, kpp, kps, None], [None, ksp, kss, ksa], [kau, None, kas, kaa]])
# ..... create linear_form
ZerVG = Constant(VG.mesh, default_scalar_type((0, 0, 0)))
Zer_AL2 = Constant(AL.mesh, default_scalar_type((0, 0, 0)))
fu = tr_S2[i] * del_u[i] * ds(2)
fp = inner(Constant(Q.mesh, 0.), del_p) * dx
fs = ZerVG[i] * del_gradu[i, j] * noder[j] * ds
fa = Zer_AL2[i] * del_lamu[i, j] * noder[j] * ds
L = form([fu, fp, fs, fa])
# ....................................... Create BCs
Bottom_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Bottom)
# ............................
clamp = Constant(V.mesh, [0., 0., 0.])
clamp_dofs = locate_dofs_topological(V, mesh.topology.dim - 1, Bottom_facets)
bc0 = dirichletbc(clamp, clamp_dofs, V)
# .........Ground BC
ground = Constant(Q.mesh, 0.0)
ground_dofs = locate_dofs_topological(Q, mesh.topology.dim - 1, Bottom_facets)
bc1 = dirichletbc(ground, ground_dofs, Q)
# ....
bcs = [bc0, bc1]
# ------------------------------  Assemble the nested system
A = assemble_matrix_nest(a, bcs=bcs)
A.assemble()
b = assemble_vector_nest(L)
apply_lifting_nest(b, a, bcs=bcs)
spaces = extract_function_spaces(L)
bcs0 = bcs_by_block(spaces, bcs)
set_bc_nest(b, bcs0)
# ------------------------------  Create a block preconditioner
paa = form(lamu[j, k] * del_lamu[j, k] * dx)
P33 = assemble_matrix(paa, [])
P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None, None, None], [None, A.getNestSubMatrix(1, 1), None, None], [None, None, A.getNestSubMatrix(2, 2), None], [None, None, None, P33]])
P.assemble()
# ------------------------------  Create Subspace solver
ksp = PETSc.KSP().create(MPI.COMM_WORLD)
ksp.setOperators(A, P)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-12, atol= 1e-12)
ksp.getPC().setType("fieldsplit")

nested_IS = P.getNestISs()
ksp.getPC().setFieldSplitIS(("u", nested_IS[0][0]),
("p", nested_IS[0][1]),
("lamu", nested_IS[0][3]))

# Set the preconditioners for each block
ksp_u, ksp_p, ksp_gradu, ksp_lamu = ksp.getPC().getFieldSplitSubKSP()

ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")

ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")

ksp_lamu.setType("preonly")
ksp_lamu.getPC().setType("gamg")
# Monitor the convergence of the KSP
ksp.setFromOptions()
###  SOLVING
Displ, potent, DisGrad, LaMA = Function(V), Function(Q), Function(VG), Function(AL)
Result = PETSc.Vec().createNest([Displ.vector, potent.vector, DisGrad.vector, LaMA.vector ])
ksp.solve(b, Result)
print(ksp.getConvergedReason())
Displ.x.scatter_forward()
potent.x.scatter_forward()