Hello everyone! When I solve the frequency response of acantilever beam, I found that the solution obtained from direct solver and CG solver were different. Even I set the rtol to 1e-20, it still didin’t work. But the same code worked when set the frequency of force to 0. Besids, when I use GMRES solver, it works for freq=50.
The follows is my MWC:
import basix
import dolfinx
import dolfinx.fem.petsc
from mpi4py import MPI
import numpy as np
import ufl
from petsc4py import PETSc
def epsilon(u):
"""
Strain function
Args:
u (dolfinx.fem.Function): Displacement function
Returns:
Strain tensor function
**Mathematical formulation:**
The strain tensor is calculated using the formula:
.. math::
\\epsilon(u) = \\frac{1}{2}\\left(\\nabla u + (\\nabla u )^T\\right)
"""
return ufl.sym(ufl.grad(u))
def sigma(u, lmbda, mu): # 3D or plane strain
"""
Cauchy stress tensor function
Args:
u (dolfinx.fem.Function): displacement function
lmbda (float | dolfinx.fem.Constant): Lamé’s elasticity parameter
mu (float | dolfinx.fem.Constant): Lamé’s elasticity parameter
Returns:
Cauchy stress tensor function
**Mathematical formulation:**
The Cauchy stress tensor (\\sigma) is calculated using the Lamé parameters and the strain tensor (\\epsilon):
.. math::
\\sigma(u) =\\lambda(\\nabla \\cdot u)I + \\mu(\\nabla u + (\\nabla u)^T)
"""
return lmbda * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
def main():
# parameters
pdic = {
"lC": 8, "wC": 2, "hC": 2, # the size of the model
'Nx':16, 'Ny': 4, 'Nz': 4, # the number of the coarse grids
"E":210e9,
"nu":0.3,
"rho":7.8e3,
"force": (1.0e5)/3, "freq": 50 # the force and frequency of the dynamic analysis
}
mu = (pdic['E'] / (2 * (1 + pdic['nu'])))
lmbda = pdic['nu'] * pdic['E'] / ((1 - 2 * pdic['nu']) * (1 + pdic['nu']))
# Create mesh-------------------------------------------------------------------------------------------------------
mesh = dolfinx.mesh.create_box(comm=MPI.COMM_WORLD, points=[[0, 0, 0], [pdic['lC'], pdic['wC'], pdic['hC']]],
n=[pdic['Nx'], pdic['Ny'], pdic['Nz']],
cell_type=dolfinx.mesh.CellType.hexahedron)
# Create the function space-----------------------------------------------------------------------------------------
c_el = basix.ufl.element(family="Lagrange", cell=mesh.topology.cell_name(),
degree=1, shape=(mesh.topology.dim,))
V = dolfinx.fem.functionspace(mesh, c_el)
# Define the boundary conditions------------------------------------------------------------------------------------
bc = dolfinx.fem.dirichletbc(np.array([0, 0, 0], dtype=dolfinx.default_scalar_type),
dolfinx.fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0)),
V)
fixed_dofs = bc._cpp_object.dof_indices()[0]
# Define the force--------------------------------------------------------------------------------------------------
p1 = lambda x: np.isclose(x[0], pdic['lC']) & np.isclose(x[1], 0) & np.isclose(x[2], 0)
p2 = lambda x: np.isclose(x[0], pdic['lC']) & np.isclose(x[1], pdic['wC'] / 2) & np.isclose(x[2], 0)
p3 = lambda x: np.isclose(x[0], pdic['lC']) & np.isclose(x[1], pdic['wC']) & np.isclose(x[2], 0)
location = lambda x: p1(x) | p2(x) | p3(x)
f_loc = dolfinx.mesh.locate_entities(mesh, 0, location)
marker = dolfinx.mesh.meshtags(mesh, 0, f_loc, np.full(len(f_loc), 1, dtype=np.int32))
f = np.zeros(mesh.geometry.x.size)
f[marker.find(1)*3+2] = -pdic['force'] # force in z-axis
# Create stiffness and mass matrices--------------------------------------------------------------------------------
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
ak = ufl.inner(sigma(u, lmbda, mu), epsilon(v)) * ufl.dx
k_form = dolfinx.fem.form(ak)
K = dolfinx.fem.petsc.assemble_matrix(k_form)
K.assemble()
am = pdic['rho'] * ufl.dot(v, u) * ufl.dx
m_form = dolfinx.fem.form(am)
M = dolfinx.fem.petsc.assemble_matrix(m_form)
M.assemble()
# dynamic stiffness matrix
w = 2 * np.pi * pdic["freq"]
W = -w ** 2 * M + K
# Apply boundary conditions using lifting method--------------------------------------------------------------------
b = PETSc.Vec().createWithArray(f)
x_ref = PETSc.Vec().createWithArray(np.zeros_like(f))
x = PETSc.Vec().createWithArray(np.zeros_like(f))
W_ref = W.copy()
# Apply lifting on rhs
g = W.createVecRight()
g.array[:] = 0
g.setValues(fixed_dofs, np.zeros_like(fixed_dofs))
g.assemble()
W_ref.multAdd(-g, b, b)
# set rhs boundary conditions
b.setValues(fixed_dofs, np.zeros_like(fixed_dofs))
# Apply lifting on matrix
W_ref.zeroRowsColumns(fixed_dofs, diag=1)
# Solve the system using Direct solver------------------------------------------------------------------------------
ksp0 = PETSc.KSP().create(comm=MPI.COMM_WORLD)
ksp0.setType(PETSc.KSP.Type.PREONLY)
ksp0.getPC().setType(PETSc.PC.Type.LU)
ksp0.setOperators(W_ref)
ksp0.setUp()
ksp0.solve(b, x_ref)
# Save the solution using cg--------------------------------------------------------------------------------------------
ksp1 = PETSc.KSP().create(comm=MPI.COMM_WORLD)
ksp1.setType(PETSc.KSP.Type.CG)
ksp1.getPC().setType(PETSc.PC.Type.ILU)
ksp1.setOperators(W_ref)
ksp1.setUp()
ksp1.solve(b, x)
print(x.array)
print(x_ref.array)
print(np.allclose(x.array, x_ref.array))
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "solution.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
funx = dolfinx.fem.Function(V)
funx.x.petsc_vec.array[:] = x_ref.array
xdmf.write_function(funx, 0.0)
if __name__ == '__main__':
main()
When pdic['freq']=0, output:
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 1.07782458e-05
-1.97558494e-08 -6.00293039e-05]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 1.07782458e-05
-1.97558494e-08 -6.00293039e-05]
True
When pdic['freq']=50, output:
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 1.76022382e-05
-1.02110686e-05 -6.96968041e-05]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... -2.13367591e-06
-1.92054057e-08 1.89030640e-05]
False
I don’t know if it’s a problem with my own code or the CG solver. Thank you very much for your help!
I will study this book.