Matrix-free for deformable solid #1
Comments
|
The first question is for @tzanio because there are private repos for internal/reporting purposes. Doesn't that dim-4 tensor have some symmetries? The crossover point may be different if you have very heavy coefficients, but the matrix-free representation is asymptotically better as the approximation order increases (and in practice the crossover tends to be pretty low, like p=2 or p=3). |
|
Here is a "back of the envelope" computation of the memory cost for matrix-free and sparse matrix (in CSR format) elasticity operator representations (produced by the python script below): Here, Unless I made a mistake, it turns out that only for Here is the python script: from numpy import *
size_int=4
size_float=8
def set_params_3d():
dim=3
p1=arange(1,3) # polynomial orders, q=p+1
p2=arange(1,17) # polynomial orders, q=p+2
p=append(p1,p2)
q=append(p1+1,p2+2) # quadrature points in 1d
qf=3**4 # number of floating point values per quadrature point
one=ones(p.shape)
nd = 3e6 # approximate (scalar) number of dofs
E=floor(nd**(1./3)/p)**3 # number of elements based on approx. num. of dofs
n=(E**(1./3)*p+1)**3 # number of scalar dofs on a E^(1/3) x E^(1/3) x E^(1/3) mesh
return dim,p1,p,q,qf,E,n
# table formatting:
fmt='%3d | %2d |%9d |%9d |%7.2f |%6.2f |%8.2f |%9.1f'
row1=' | | num. of | num. of | T-vec | ratio | B/dof | B/dof'
row2=' p | q | elements | dofs | MiB | #E/#T | MF | CSR'
sep= '----+----+----------+----------+--------+-------+---------+----------'
big_sep='='*70
def print_table(t):
print row1
print row2
print sep
for i in range(len(t)):
if i==len(p1):
print sep
print fmt%tuple(t[i])
# ----------------------------------------------------------
# HEXES ELASTICITY
# ----------------------------------------------------------
def elast_hex_mf_bytes(p,q,qf,E):
# Element-to-dof array + quadrature point data
return (E*dim*(p+1)**3*size_int+(E*q**3)*qf*size_float)
def elast_hex_csr_bytes(p,E):
e=E**(1./3)
# csr nnz on a periodic (e x e x e) mesh
# num vert dofs x num vert dof connections:
nnz=dim*(e+1)**3 * dim*(2*p+1)**3
# add num edge dofs x num edge dof connections:
nnz=nnz + dim*(p-1)*(dim*e*(e+1)**2) * dim*(2*p+1)**2*(p+1)
# add num face dofs x num face dof connections:
nnz=nnz + dim*(p-1)**2*(dim*e**2*(e+1)) * dim*(2*p+1)*(p+1)**2
# add num element dofs x num element dof connections:
nnz=nnz + dim*(p-1)**3*(e**3) * dim*(p+1)**3
# 1 int + 1 float for every nonzero entry + (mat_size+1) ints
return nnz*(size_int+size_float)+(dim*n+1)*size_int
dim,p1,p,q,qf,E,n=set_params_3d()
print big_sep
print 'elasticity operator action, dimension:', dim
c5=dim*n*size_float/2.**20 # T-vector size in MiB
c6=E*(p+1)**3/(1.*n) # E-vector size / T-vector size
c7=elast_hex_mf_bytes(p,q,qf,E)/(1.*dim*n) # mf bytes per dof
c8=elast_hex_csr_bytes(p,E)/(1.*dim*n) # csr bytes per dof
t=array([p,q,E,dim*n,c5,c6,c7,c8]).transpose()
print
print_table(t)
print big_sep |
|
Hi Bruno, Regarding your first question, we will discuss internally if everyone is comfortable with making the milestone reports public and get back to you. We do try to highlight the major developments on our News page, http://ceed.exascaleproject.org/news/, but I agree that the technical content in the reports can be interesting and valuable to an external audience. Cheers, |
|
@v-dobev Thanks for the estimation and for the script! This is actually much better than what we expected. They said that they will try to use matrix-free in their application. We will keep you posted. Thanks again. @tzanio Thanks. |
|
This figure plots cost estimates with a common symmetry symmetry. |
|
Here is a more realistic version of the CSR storage computation where the sparsity pattern for one scalar components is reused for all 3 x 3 sparse blocks: Note that still the cheapest CSR case (p=1) uses 675.9 bytes/dof while a moderate order (p=5) MF needs less: 586.95 bytes/dof. If you go to order p=12, MF is about 2x cheaper in terms of memory. Here is the updated python script: from numpy import *
size_int=4
size_float=8
def set_params_3d():
dim=3
p1=arange(1,3) # polynomial orders, q=p+1
p2=arange(1,17) # polynomial orders, q=p+2
p=append(p1,p2)
q=append(p1+1,p2+2) # quadrature points in 1d
qf=3**4 # number of floating point values per quadrature point
one=ones(p.shape)
nd = 3e6 # approximate (scalar) number of dofs
E=floor(nd**(1./3)/p)**3 # number of elements based on approx. num. of dofs
n=(E**(1./3)*p+1)**3 # number of scalar dofs on a E^(1/3) x E^(1/3) x E^(1/3) mesh
return dim,p1,p,q,qf,E,n
# table formatting:
fmt='%3d | %2d |%9d |%9d |%7.2f |%6.2f |%8.2f |%9.1f'
row1=' | | num. of | num. of | T-vec | ratio | B/dof | B/dof'
row2=' p | q | elements | dofs | MiB | #E/#T | MF | CSR'
sep= '----+----+----------+----------+--------+-------+---------+----------'
big_sep='='*70
def print_table(t):
print row1
print row2
print sep
for i in range(len(t)):
if i==len(p1):
print sep
print fmt%tuple(t[i])
# ----------------------------------------------------------
# HEXES ELASTICITY
# ----------------------------------------------------------
def elast_hex_mf_bytes(p,q,qf,E):
# Element-to-dof array + quadrature point data
return (E*dim*(p+1)**3*size_int+(E*q**3)*qf*size_float)
def elast_hex_csr_bytes(p,E,n):
e=E**(1./3)
#
# csr nnz on a (e x e x e) hex mesh
# (one scalar component, full local matrices)
#
# dofs on vertices with one adjacent element:
nnz=2**3 * (p+1)**3
# dofs on vertices with two adjacent elements:
nnz=nnz + 12*(e-1) * (2*p+1)*(p+1)**2
# dofs on vertices with four adjacent elements:
nnz=nnz + 6*(e-1)**2 * (2*p+1)**2*(p+1)
# dofs on vertices with eight adjacent elements:
nnz=nnz + (e-1)**3 * (2*p+1)**3
#
# dofs on edges with one adjacent element:
nnz=nnz + (p-1)*12*e * (p+1)**3
# dofs on edges with two adjacent elements:
nnz=nnz + (p-1)*6*2*e*(e-1) * (2*p+1)*(p+1)**2
# dofs on edges with four adjacent elements:
nnz=nnz + (p-1)*3*e*(e-1)**2 * (2*p+1)**2*(p+1)
#
# dofs on faces with one adjacent element:
nnz=nnz + (p-1)**2*6*e**2 * (p+1)**3
# dofs on faces with two adjacent elements:
nnz=nnz + (p-1)**2*3*e**2*(e-1) * (2*p+1)*(p+1)**2
#
# dofs inside elements:
nnz=nnz + (p-1)**3*e**3 * (p+1)**3
#
# 1 int + 1 float for every nonzero entry + (mat_size+1) ints
# return 9*nnz*(size_int+size_float)+(3*n+1)*size_int
#
# reuse the scalar sparsity pattern for all 9 blocks:
return 9*nnz*size_float + nnz*size_int + (n+1)*size_int
#
# ignoring any sparsity pattern storage
# return 9*nnz*size_float
dim,p1,p,q,qf,E,n=set_params_3d()
print big_sep
print 'elasticity operator action, dimension:', dim
c5=dim*n*size_float/2.**20 # T-vector size in MiB
c6=E*(p+1)**3/(1.*n) # E-vector size / T-vector size
c7=elast_hex_mf_bytes(p,q,qf,E)/(1.*dim*n) # mf bytes per dof
c8=elast_hex_csr_bytes(p,E,n)/(1.*dim*n) # csr bytes per dof
t=array([p,q,E,dim*n,c5,c6,c7,c8]).transpose()
print
print_table(t)
print big_sep |
I just wanted to echo @Rombur's thanks to both @v-dobrev and @jedbrown. This led to an interesting discussion on our side, and as @Rombur said we'll be investigating this more and will report back when we have some results to share. |
|
You are welcome! Glad this is helpful. |
|
The CEED reports are now available here: |
|
@Rombur @jppelteret With regard to solid mechanics and efficient storage for the matrix-free Jacobian, please check out the discussion in this section of the libCEED docs, including the note on storage complexity below. https://libceed.readthedocs.io/en/latest/examples/solids/#id3 |
|
We finally investigated different caching strategies in https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6336. For the example considered there, it was best to just cache a second-rank tensor while caching the full rank-4 tensor was slowest. |
|
Nice paper @masterleinad. It looks like Table 1 Algorithm 3 isn't exploiting symmetry in the fourth order elasticity tensor C, would make it 21 instead of 36 independent entries (in 3D). Our currently-implemented formulation is similar to the discussion at the end of your appending, but (I think) exploiting more symmetry. Fusing the metric transformation to the current configuration for Jacobian application (as in your Algorithm 2) is also clearly a good choice. It'll be interesting to make quantitative comparisons. Is your code available? |
|
The code can be found at https://github.com/davydden/large-strain-matrix-free/. |
|
Thanks @jedbrown for the link to the CEED documentation. That example and writeup looks great! And thanks once more to you and @v-dobrev for the initial discussions on this topic. I was going to pass along the link to the paper, but @masterleinad seems to have beaten me to it :-)
This is the current implementation in of symmetric tensors deal.II (they lack major symmetry). Implementing a new tensor class was a little beyond the scope of the work, so we didn't consider it at all.
Great that you got that covered already :-) |
|
@jedbrown I'd be very happy in looking into your implementations as well; it's been an interesting read through your materials and the general overview - let me know when you start looking into our algorithms so I can dig deeper into yours and make cross-comparisons. In any case, a discussion from both teams might provide further insight into what we are doing, so 👍 |
|
Indeed! We'd like to refine a couple things in our implementation and set up some good comparisons. We've been doing experiments with load increments (which have lots of opportunities for nonlinear multilevel methods), but make direct comparison a bit tricky. We'll have a look at your code and see if we can set up a fair comparison. One thing we noticed that made a significant impact for the stored-nothing Jacobian method was to replace Our implementation will be ready for GPU comparisons once we can compute the diagonal of the Jacobian efficiently there. Overall, I think the approaches that your group and ours are using should have significant impact on industrial practice, where assembled matrices (often with quadratic elements) are the norm. |
First, I have a general comment. Would it be possible to post here or on the website, the milestone reports that are produced as part of ECP. Right now it is hard for people not affiliated with ECP to follow the developments made by CEED. While someone could, in theory, stalk all the repositories that are affiliated with CEED, it would be a lot easier to just read the milestone reports.
Second, I have a question about the limitation of matrix-free methods. I have some colleagues simulating the deformation of solid materials with complex constitutive laws. In this case, computing the material properties on the fly is extremely costly due to the complexity of these laws that sometimes require the inversion of a small linear system. Yet precomputing these properties is not satisfactory either because we can end up storing a fourth-order tensor (about 80 coefficients) per quadrature points. This means that the matrix-free implementation requires more memory than the matrix-based implementation. We could of course try to store some elements of the constitutive laws and compute the rest on the fly but it is not obvious how to strike the right balance. Is there anyone in CEED who is facing a similar problem?
cc: @bangerth, @tjhei, @kronbichler, @masterleinad, @davydden, @jppelteret
The text was updated successfully, but these errors were encountered: