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

arash_kazemi · May 17, 2024, 9:50pm

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))
# ......................................... Reading mesh file
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
# ......................................... Define Loading
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])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bc)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

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)

francesco-ballarin · May 18, 2024, 5:22am

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

Apart from this, it’s impossible for anyone to help you any further, because we don’t have access to the meshes and thus we can’t run your code.

arash_kazemi · May 18, 2024, 7:41pm

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


import gmsh
import sys

##################
gmsh.initialize()
gmsh.model.add("Model")
# =====================================    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)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 12, (LN - 3), -(LN - 8), -(LN - 4)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 11, (LN - 2), -(LN - 7), -(LN - 3)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 10, (LN - 1), -(LN - 6), -(LN - 2)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 9, (LN - 4), -(LN - 5), -(LN - 1)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 8, LN - 7, LN - 6, LN - 5], SN)
gmsh.model.geo.addPlaneSurface([SN], 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)
gmsh.model.geo.addVolume([VN], VN)
Vol_list.append(VN)
VN += 1
gmsh.model.geo.addPoint(Points[i][0], 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], 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)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list2.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 15, LN - 6, -(LN - 3), -(LN - 7)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list3.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 14, LN - 5, -(LN - 2), -(LN - 6)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list4.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 13, LN - 8, -(LN - 1), -(LN - 5)], SN)
gmsh.model.geo.addPlaneSurface([SN], 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], [])

gmsh.model.geo.addPoint(Left_down_corner[0], Left_down_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Right_down_corner[0], Right_down_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Right_up_corner[0], Right_up_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Left_up_corner[0], Left_up_corner[1], 0, lc, PN)
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)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list2.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 11, LN - 6, -(LN - 3), -(LN - 7)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list3.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 10, LN - 5, -(LN - 2), -(LN - 6)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list4.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 9, LN - 8, -(LN - 1), -(LN - 5)], SN)
gmsh.model.geo.addPlaneSurface([SN], 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], [])
gmsh.model.geo.addPoint(Left_down_corner[0], Left_down_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Right_down_corner[0], Right_down_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Right_up_corner[0], Right_up_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Left_up_corner[0], Left_up_corner[1], 0, 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)
# ################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)
gmsh.model.geo.addPlaneSurface([SN], 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)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list2.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 27, LN - 19, (LN - 11), (LN - 7)], SN)
gmsh.model.geo.addPlaneSurface([SN], 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)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list2.append(SN - 1)

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

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

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

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

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

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

gmsh.model.geo.addSurfaceLoop(Sur_list5, VN)
gmsh.model.geo.addVolume([VN], VN)
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)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list6.append(SN - 1)

gmsh.model.geo.addSurfaceLoop(Sur_list6, VN)
gmsh.model.geo.addVolume([VN], VN)
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.model.addPhysicalGroup(3, Vol_list, name="Pyramid")
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()

arash_kazemi · May 20, 2024, 8:32pm

@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

francesco-ballarin · May 21, 2024, 5:40am

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?

arash_kazemi · May 21, 2024, 5:35pm

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,

francesco-ballarin · May 21, 2024, 5:59pm

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.

arash_kazemi · May 21, 2024, 6:49pm

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

arash_kazemi · May 31, 2024, 6:40pm

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,

dokken · May 31, 2024, 7:47pm

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?

arash_kazemi · May 31, 2024, 10:26pm

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?

arash_kazemi · June 3, 2024, 6:21pm

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

Thanks

dokken · June 3, 2024, 7:36pm

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.

arash_kazemi · June 3, 2024, 8:30pm

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

arash_kazemi · June 19, 2024, 3:34am

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
# ......................................... Reading mesh file
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)
# ......................................... Define Loading
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])


def RStrianGrad(gradu):
    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
kps = -(flexo[i, j, k, l] * RStrianGrad(gradu)[j, k, l]) * (-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
kss = (SGElas(la, non_local, mu)[i, j, k, l, m, n] * RStrianGrad(gradu)[l, m, n]) * 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()
DisGrad.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

dokken · June 19, 2024, 6:35am

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

arash_kazemi · June 19, 2024, 5:39pm

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

dokken · June 19, 2024, 7:42pm

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.

arash_kazemi · June 23, 2024, 7:02pm

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
# ......................................... Reading mesh file
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)
# ......................................... Define Loading
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
def RStrianGrad(gradu):
    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
kps = -(flexo[i, j, k, l] * RStrianGrad(gradu)[j, k, l]) * (-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
kss = (SGElas(la, non_local, mu)[i, j, k, l, m, n] * RStrianGrad(gradu)[l, m, n]) * 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")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)

nested_IS = P.getNestISs()
ksp.getPC().setFieldSplitIS(("u", nested_IS[0][0]),
                            ("p", nested_IS[0][1]),
                            ("gradu", nested_IS[0][2]),
                            ("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_gradu.setType("preonly")
ksp_gradu.getPC().setType("gamg")

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()
DisGrad.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)

Topic		Replies	Views
Dynamic Elasticity Not Converging dolfinx	0	166	March 1, 2024
Problems with convergence of the elasticity equation General	1	277	September 29, 2023
Discrepancy Between Direct and Iterative Solvers General	5	111	June 28, 2024
Hyperelasticity example issues General	9	209	May 3, 2024
3D elasticity model with a big mesh Errors	0	570	September 29, 2022

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

Related Topics