hi,I tried to repeat this plastic numerical example, 2.077-Fenicsx-Plasticity/Plasticity_Tests/2d_test_Plasticity_New_Frac_Snes.ipynb at main · jorgenin/2.077-Fenicsx-Plasticity · GitHub but there was a problem. I simplified the problem into a simple elastic problem, and the following error message occurred. May I ask why this is the case? The solver reference is the same as that used in the numerical example above.
x = int(1 / 0.1)
ny = int(1 / 0.1)
domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, cell_type=dolfinx.mesh.CellType.triangle)MPI.COMM_WORLD, 0, 2)
tdim =2
fdim =1
V = fem.functionspace(domain,("CG",1,(2,)))
u_new = fem.Function(V,name = "Displacement_for_updating")
du_limit = fem.Function(V, name="Displacement_limit")
du_max = fem.Function(V, name="Displacement_max")
u_load = fem.Function(V, name="Displacement_load")
v = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
du_limit.x.array[:] = -np.inf
du_max.x.array[:] = np.inf
n = ufl.FacetNormal(domain)
loading = fem.Constant(domain,ScalarType(0.0))
def fix(x):
return np.isclose(x[1],0.0)
def load(x):
return np.isclose(x[1],1.0)
def F_ext(v):
return -loading * ufl.inner(n, v) * ds
BC = []
Uimp = fem.Constant(domain,1.0)
dofs_fix = dolfinx.fem.locate_dofs_geometrical(V,fix)
entities = dolfinx.mesh.locate_entities_boundary(domain,fdim,load)
dofs_load = dolfinx.fem.locate_dofs_topological(V.sub(1),fdim,entities)
load_tag = mesh.meshtags(domain,fdim,entities,np.full(len(entities),2))
BC_fix = dolfinx.fem.dirichletbc(dolfinx.default_scalar_type((0.0,0.0)),dofs_fix,V)
BC_load = dolfinx.fem.dirichletbc(Uimp,dofs_load,V.sub(1))
BC = [BC_fix]
u_load.x.array[:] = Uimp.value[()]
E = fem.Constant(domain, 7e3)
nu = fem.Constant(domain, 0.1)
mu = E / (2.0 * (1.0 + nu))
lamda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
def eps(u):
return ufl.sym(ufl.grad(u))
def sigma(u):
return lamda * ufl.tr(eps(u)) * ufl.Identity(2) + 2 * mu * eps(u)
dx = ufl.Measure("dx",domain=domain)
ds = ufl.Measure("ds",domain = domain,subdomain_data=load_tag)
E_u = ufl.inner(sigma(u_new),eps(v)) * dx
E_u_jacbian = ufl.derivative(E_u,u_new,u_)
problem = SNESSolver(E_u, u_new,J_form = E_u_jacbian, bcs = BC, petsc_options=petsc_options_SNES)
savedir = "results1/"
if os.path.isdir(savedir):
shutil.rmtree(savedir)
os.makedirs(savedir + "Displacement", exist_ok=True)
Displacement_files = dolfinx.io.VTKFile(domain.comm, savedir + "Displacement/Displacement.vtk", "w")
Nincr = 50
loading = np.linspace(0, 0.01,Nincr + 1)
Niter_max = 100
tol = 1e-4
tol_plastic = 1e-4
iterations = [0]
Force = [0]
for iter, t in enumerate(loading[1:]):
print("Load step", iter + 1)
Uimp.value[()] = t
problem.solve()
iterations.append(iter)
Force.append(fem.assemble_scalar(fem.form(sigma(u_new)[1,1] *ds(2))))
Displacement_files.write_function(u_new, iter)
plt.figure()
plt.plot(loading, Force)
plt.xlabel("Imposed displacement")
plt.ylabel("Reaction force")
plt.show()
the snes solver
from dolfinx import mesh, fem, io
import dolfinx.fem.petsc
from petsc4py import PETSc
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import os, shutil
from petsc4py.PETSc import ScalarType
from dolfinx.cpp.log import LogLevel, log
petsc_options_SNES = {
"snes_type": "vinewtonrsls",
"snes_linesearch_type": "basic",
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"snes_atol": 1.0e-08,
"snes_rtol": 1.0e-09,
"snes_stol": 0.0,
"snes_max_it": 50,
"snes_monitor": "",
# "snes_monitor_cancel": "",
}
class SNESSolver:
"""
Problem class for elasticity, compatible with PETSC.SNES solvers.
"""
def __init__(
self,
F_form: ufl.Form,
u: fem.Function,
bcs=[],
J_form: ufl.Form = None,
bounds=None,
petsc_options={},
form_compiler_parameters={},
jit_parameters={},
monitor=None,
prefix=None,
):
self.u = u
self.bcs = bcs
self.bounds = bounds
# Give PETSc solver options a unique prefix
if prefix is None:
prefix = "snes_{}".format(str(id(self))[0:4])
self.prefix = prefix
if self.bounds is not None:
self.lb = bounds[0]
self.ub = bounds[1]
V = self.u.function_space
self.comm = V.mesh.comm
self.F_form = fem.form(F_form)
if J_form is None:
J_form = ufl.derivative(F_form, self.u, ufl.TrialFunction(V))
self.J_form = fem.form(J_form)
self.petsc_options = petsc_options
self.monitor = monitor
self.solver = self.solver_setup()
def set_petsc_options(self, debug=False):
# Set PETSc options
opts = PETSc.Options()
opts.prefixPush(self.prefix)
if debug is True:
ColorPrint.print_info(self.petsc_options)
for k, v in self.petsc_options.items():
opts[k] = v
opts.prefixPop()
def solver_setup(self):
# Create nonlinear solver
snes = PETSc.SNES().create(self.comm)
# Set options
snes.setOptionsPrefix(self.prefix)
self.set_petsc_options()
snes.setFromOptions()
self.b = fem.petsc.create_vector(self.F_form)
self.a = fem.petsc.create_matrix(self.J_form)
snes.setFunction(self.F, self.b)
snes.setJacobian(self.J, self.a)
# We set the bound (Note: they are passed as reference and not as values)
if self.monitor is not None:
snes.setMonitor(self.monitor)
if self.bounds is not None:
snes.setVariableBounds(self.lb.vector, self.ub.vector)
return snes
def F(self, snes: PETSc.SNES, x: PETSc.Vec, b: PETSc.Vec):
"""Assemble the residual F into the vector b.
Parameters
==========
snes: the snes object
x: Vector containing the latest solution.
b: Vector to assemble the residual into.
"""
# We need to assign the vector to the function
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
x.copy(self.u.x.petsc_vec)
self.u.x.scatter_forward()
# Zero the residual vector
b.array[:] = 0
fem.petsc.assemble_vector(b, self.F_form)
# this is a nasty workaround to include the force term with the bug https://github.com/FEniCS/dolfinx/issues/2664
force_form = fem.form(-F_ext(v))
b_ds = fem.petsc.create_vector(force_form)
fem.petsc.assemble_vector(b_ds,force_form)
b.array[:] += b_ds.array
# Apply boundary conditions
fem.petsc.apply_lifting(b, [self.J_form], [self.bcs], [x], -1.0)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(b, self.bcs, x, -1.0)
def J(self, snes, x: PETSc.Vec, A: PETSc.Mat, P: PETSc.Mat):
"""Assemble the Jacobian matrix.
Parameters
==========
x: Vector containing the latest solution.
A: Matrix to assemble the Jacobian into.
"""
A.zeroEntries()
fem.petsc.assemble_matrix(A, self.J_form, self.bcs)
A.assemble()
def solve(self):
log(LogLevel.INFO, f"Solving {self.prefix}")
try:
self.solver.solve(None, self.u.x.petsc_vec)
self.u.x.scatter_forward()
return (self.solver.getIterationNumber(), self.solver.getConvergedReason())
except Warning:
log(
LogLevel.WARNING,
f"WARNING: {self.prefix} solver failed to converge, what's next?",
)
raise RuntimeError(f"{self.prefix} solvers did not converge")
type or paste code here
the problem would be like that
b_ds = fem.petsc.create_vector(force_form)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 112, in create_vector
dofmap = L.function_spaces[0].dofmap
^^^^^^^^^^^^^^^^^
AttributeError: 'list' object has no attribute 'function_spaces'
type or paste code here
