I’m solving a 3D finite element problem where the number of degrees of freedom of the mesh exceeds 2^{16}, which is the representation range of int32. The code reports the following error:
u_f_vec.setValues(imap, Sele_disp[i], addv=PETSc.InsertMode.INSERT_VALUES)
File "PETSc/Vec.pyx", line 1025, in petsc4py.PETSc.Vec.setValues
File "PETSc/petscvec.pxi", line 332, in petsc4py.PETSc.vecsetvalues
File "include/arraynpy.pxi", line 137, in petsc4py.PETSc.iarray_i
File "include/arraynpy.pxi", line 130, in petsc4py.PETSc.iarray
TypeError: Cannot cast array data from dtype('uint32') to dtype('int32') according to the rule 'safe'
I looked up the petsc4py manual and found that PETSc has two versions 64bit and 32bit (default). I guess for my problem I may need to use the 64bit version of petsc4py.I searched the forums and found that the only relevant posting known is installing the 64bit version on Docker. But I am using Ubuntu, so I would like to know how to install and configure 64bit petsc4py on Ubuntu.
Thank you for your suggestions. But the max index of dofs for my problem is exceed the range that the type np.int32 can represent. So I must use a longer type to represent the integer.
You would have to first install the specified petsc (and uninstall he other petsc). I al not sure it will be picked up correctly though, as it wasnt listed in the dependencies on the packages page.
Thanks Mr.dokken. And I still has a question to ask you for help.
I read the manual of PETSc, and it reads that:
[When Should/Can I Use The configure Option --with-64-bit-indices?(Frequently Asked Questions (FAQ) — PETSc 3.15.0 documentation)
By default the type that PETSc uses to index into arrays and keep sizes of arrays is a PetscInt defined to be a 32 bit int. If your problem:
Involves more than 2^31 - 1 unknowns (around 2 billion).
Your matrix might contain more than 2^31 - 1 nonzeros on a single process.
Then you need to use this option. Otherwise you will get strange crashes.
This option can be used when you are using either 32 bit or 64 bit pointers. You do not need to use this option if you are using 64 bit pointers unless the two conditions above hold.]
The max value of degree of freedom n for my problem is 2^{16}<n<2^{31}-1, which is less than 2 billon. But when setting the value of the Vec, the indices type must be np.int32 in the method Vec.setValues(indices, values[, addv]). The largest value of np.int32 can represent is 2^{16}, so I wouldn’t be able to set the values of Vec in my problem.
My program is serial right now, and I may consider using parallelism later. Because PETSc is known to have the ability to support massively parallel computing I wanted to develop my program based on it. With regard to the issues mentioned above, let me illustrate with a simple example. For example I want to apply a concentrated force or a distributed force along a certain edge by modifying the vector value of the right end term of the balance equation.
def statics_solver_petsc(self, f, K, fixed, petsc_options={"ksp_type": "PREONLY", "pc_type": "LU"}):
"""
Solve the static equilibrium equation.
Args:
f (np.ndarray): the force array
K (PETSc.Mat): the stiffness matrix
fixed (np.ndarray): the dofs that is fixed
petsc_options (dict):
Returns:
"""
# Create the ksp solver-----------------------------------------------------------------------------------------
ksp_type = getattr(PETSc.KSP.Type, petsc_options["ksp_type"])
pc_type = getattr(PETSc.PC.Type, petsc_options["pc_type"])
ksp = PETSc.KSP()
ksp.create(comm=self.comm)
ksp.setType(ksp_type)
ksp.getPC().setType(pc_type)
ksp.getPC().setFactorSolverType("mumps") # Set the solver package used to perform the factorization.
# Set boundary condition----------------------------------------------------------------------------------------
K.zeroRows(fixed, diag=1)
ksp.setOperators(K)
x, b = K.createVecs() # Create a displacement vector and a load vector
b.setValues(np.arange(b.getSize(), dtype=np.uint32), f)
# Solve the problem
ksp.solve(b, x)
return x
The above problem occurs in the process of modifying the Vec variable b. Because here I’m using a 32-bit unsigned integer np.uint32 to represent the degrees of freedom indices, which is outside the range of np.int32.
You should never have that many degrees of freedom on a single process. If you want to run such large problems, you have to think about how to run it in a distributed sense.
Could you please create a minimal reproducible example, where you use a built in grid (unit cube mesh?) to assemble a stiffness matrix over, and then show what you want to change (as atm f and K are not defined in your code.
Thanks for your advice. It is indeed bad to include too many degrees of freedom on a single process, and the best way is to use distributed memory and parallel computing. However I haven’t been exposed to FEniCSx and PETSc, and parallel computing for long enough to know how to utilize distributed memory and parallel computing. So the implementation was done using serial computing first. But later I went back to learn how to utilize parallel computing and distributed memory.
The following code is a minimal reproducible example.
import dolfinx
import dolfinx.fem.petsc
from petsc4py import PETSc
import mpi4py
import numpy as np
import ufl
def epsilon(u):
return ufl.sym(ufl.grad(u))
def sigma(u, lmbda, mu):
return lmbda * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
def statics_solver_petsc(f, K, fixed, petsc_options={"ksp_type": "PREONLY", "pc_type": "LU"}):
"""
Solve the static equilibrium equation.
Args:
f (np.ndarray): the force array
K (PETSc.Mat): the stiffness matrix
fixed (np.ndarray): the dofs that is fixed
petsc_options (dict):
Returns:
"""
# Create the ksp solver-----------------------------------------------------------------------------------------
ksp_type = getattr(PETSc.KSP.Type, petsc_options["ksp_type"])
pc_type = getattr(PETSc.PC.Type, petsc_options["pc_type"])
ksp = PETSc.KSP()
ksp.create(comm=mpi4py.MPI.COMM_WORLD)
ksp.setType(ksp_type)
ksp.getPC().setType(pc_type)
ksp.getPC().setFactorSolverType("mumps") # Set the solver package used to perform the factorization.
# Set boundary condition----------------------------------------------------------------------------------------
K.zeroRows(fixed, diag=1)
ksp.setOperators(K)
x, b = K.createVecs() # Create a displacement vector and a load vector
b.setValues(np.arange(b.getSize(), dtype=np.uint32), f)
# Solve the problem
ksp.solve(b, x)
return x
pdic = {
"E": 1, "nu": 0.3,
}
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=mpi4py.MPI.COMM_WORLD, points=[[0, 0, 0], [8, 2, 2]], n=[160, 40, 40],
cell_type=dolfinx.mesh.CellType.hexahedron)
# Create functionspace
element = ufl.VectorElement(family="CG", cell=mesh.ufl_cell(), degree=1)
V = dolfinx.fem.functionspace(mesh=mesh, element=element)
# dirichlet bc
left_face = dolfinx.fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0.0))
bc = dolfinx.fem.dirichletbc(dolfinx.default_scalar_type((0, 0, 0)), left_face, V)
fixed, _ = bc._cpp_object.dof_indices()
# force vector
def rd_center(x):
return np.isclose(x[0], 8.0) & np.isclose(x[1], 1.0) & np.isclose(x[2], 0.0)
rdc_node = dolfinx.fem.locate_dofs_geometrical(V, rd_center)
force = np.zeros(mesh.geometry.x.size)
force[rdc_node * 3 + 2] = -1 # the concentrate force in z-direction
# variation forms
f = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0, 0, 0)))
u = ufl.TrialFunction(V) # Define trial function
v = ufl.TestFunction(V) # Define test function
a = (ufl.inner(sigma(u, lmbda, mu), epsilon(v)) * ufl.dx) # Define bilinear form
L = ufl.dot(f, v) * ufl.dx # define linear form
a_compiled = dolfinx.fem.form(a)
L_compiled = dolfinx.fem.form(L)
# assemble_matrix with Dirichlet conditions
K = dolfinx.fem.petsc.assemble_matrix(a_compiled)
K.assemble()
u_sol = statics_solver_petsc(force, K, fixed)
It reports errors when running:
Traceback (most recent call last):
File "/mnt/d/BaiduNetdiskWorkspace/MultiscaleDynamicTOP/code/FEniCSBasedProject/3D_test/demo.py", line 85, in <module>
u_sol = statics_solver_petsc(force, K, fixed)
File "/mnt/d/BaiduNetdiskWorkspace/MultiscaleDynamicTOP/code/FEniCSBasedProject/3D_test/demo.py", line 40, in statics_solver_petsc
b.setValues(np.arange(b.getSize(), dtype=np.uint32), f)
File "PETSc/Vec.pyx", line 1025, in petsc4py.PETSc.Vec.setValues
File "PETSc/petscvec.pxi", line 332, in petsc4py.PETSc.vecsetvalues
File "include/arraynpy.pxi", line 137, in petsc4py.PETSc.iarray_i
File "include/arraynpy.pxi", line 130, in petsc4py.PETSc.iarray
TypeError: Cannot cast array data from dtype('uint32') to dtype('int32') according to the rule 'safe'
Process finished with exit code 1
As you can see from the error message, this has nothing to do with the bound, but that you are sending in unsigned integers, when petsc expects signed integers.
Simply send in:
You could use libpetsc64-real-dev instead of petsc64-dev.
Likewise for 64-bit slepc/slepc4py if you need them.
You’ll want to set PETSC_DIR=/usr/lib/petsc64 to activate use of 64-bit petsc
(likewise SLEPC_DIR=/usr/lib/slepc64)
If you want finer control of which petsc build you’re using (e.g. complex rather than real), you’d point these env variables at the version you need under /usr/lib/petscdir/ (likewise /usr/lib/slepcdir/)
You are right, Mr. dokken. I made the stupid mistake of thinking that the range of np.int32 is [-2^{16}-1,2^{16}-1] , so I use np.uint32 to get a wider range of the indices. However, the range of np.int32 is actually [-2^{31}-1,2^{31}-1], after my actual verification. So the code runs now.
