Hi all,
I’m using DOLFINx version 0.8.0 to solve a coupled system in a beam model, but I’m facing convergence issues with high aspect ratio geometries. For the case of only linear elasticity problem in a beam model, with an aspect ratio of 2 (width = 50, length = 100), the simulation converges in 9.4 seconds with 11 iterations. However, when the aspect ratio increases to 20 (length = 1000), the iterative solver runs indefinitely without converging. I tried a denser mesh, but it had no effect.
The code works fine with a direct solver, but due to problem size, I need to use an iterative solver.
Here is a simplified code for a linear elasticity problem:
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, Expression)
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, create_box, CellType, exterior_facet_indices
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
from dolfinx import la
import dolfinx
print(
f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")
start = time.time()
################# Material Properties
E = 1 * 190e9
nu = 0.3
laa = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
# ......................................... Reading mesh file
H = 50
L = 1000
mesh = create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, H, H])],[15, 15, 15], cell_type=CellType.tetrahedron)
# ......................................... Define Loading
Force = 10
Area = H * H
tr = Force / Area
######################################################
tr_top = Constant(mesh, ScalarType((0, 0, -tr)))
# ......................................... Define Elements
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 = 1e-14
def Left(x):
return np.isclose(x[0], 0, atol=tol)
# ....................................... ds
boundaries = [(1, lambda x: np.isclose(x[0], 0.0, atol=tol)),
(2, lambda x: np.isclose(x[0], L, atol=tol)),
(3, lambda x: np.logical_and(np.isclose(x[0], 0, atol=tol), np.isclose(x[2], 0, 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_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Left)
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]
# ......................................... Define tensor
delta = Identity(len(u))
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(laa * 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(laa, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])
# ......................................... Define Variational Formulation
a = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
L = tr_top[i] * del_u[i] * ds(2)
#.................................................SOLVER
########## Iterative SOLVER
if mesh.comm.rank == 0:
print("Assembling...")
bilinear_form = form(a)
linear_form = form(L)
A = assemble_matrix(bilinear_form, bcs=bc)
A.assemble()
##############
b = create_vector(linear_form)
with b.localForm() as b_loc:
b_loc.set(0)
assemble_vector(b, linear_form)
apply_lifting(b, [bilinear_form], [bc])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # if run in parrallel
set_bc(b, bc)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
##########################################
Result = Function(Space)
###################################################
def monitor(ksp, it, rnorm):
if mesh.comm.rank == 0:
print(f"iteration: {it} rnorm: {rnorm:.3e}")
solver = PETSc.KSP().create(mesh.comm) # type: ignore
# # Set solver options
opts = PETSc.Options() # type: ignore
# #######
opts["ksp_type"] = "gmres"
opts["ksp_rtol"] = 1e-4
opts["ksp_max_it"] = 70000
opts["pc_type"] = "hypre"
opts["pc_hypre_type"] = "boomeramg"
#############################
solver.setFromOptions()
solver.setOperators(A)
solver.setMonitor(monitor)
solver.solve(b, Result.vector)
if mesh.comm.rank == 0:
print(solver.getConvergedReason())
###################
Result.x.scatter_forward()
# ......................................... Export Results
with VTXWriter(mesh.comm, ("Beam/DIS.bp"), [Result]) as vtx:
vtx.write(0.0)
end = time.time()
if mesh.comm.rank == 0:
print('Time', end - start)
Thanks