Solver for 2D elasticity problem

Linear Algebra
1

Hello,

The same code I was using in Fenics 2017 for solving a linear elasticity, using a LU solver, is much slower using Fenics 2019. You can see the code below. If you run it using a 102x51 mesh, then the computational times are similar. However, with a mesh 202x101, the current running time is around 25s while it was around 8s with Fenics 2017. With a 302x151 mesh, I get a UMFPACK error: out of memory, but with Fenics 2017 it was taking around 20s to solve.

Any idea what could be the origin of this issue, and what I could try to improve the computational time?

import time    
from dolfin import *    
import numpy as np

# ----------------------------------------------------------------------
def _main():

    # Define mesh
    Nx,Ny = [202,101]   # number of divisions in the x- and y-directions 
    lx,ly = [2.0,1.0]    
    XX,YY = np.meshgrid(np.linspace(0.0,lx,Nx+1),np.linspace(0.0,ly,Ny+1))   
    mesh  = RectangleMesh(Point(0.0,0.0),Point(lx,ly),Nx,Ny,'crossed') # each square of the grid is divided in four triangles
    Vvec = VectorFunctionSpace(mesh, 'CG', 1)         

    # Define boundary conditions
    class DirBd(SubDomain):
        def inside(self, x, on_boundary):
            return near(x[0],.0)      
    dirBd = DirBd()
    boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
    boundaries.set_all(0)    
    dirBd.mark(boundaries, 1)
    ds = Measure("ds")(subdomain_data=boundaries)  
    bcd  = [DirichletBC(Vvec, (0.0,0.0), boundaries, 1)]     

    # Define level set function phi    
    phi_mat = -np.cos(8.0/lx*pi*XX) * np.cos(4.0*pi*YY) - 0.4\
         + np.maximum(200.0*(0.01-XX**2-(YY-ly/2)**2),.0)\
         + np.maximum(100.0*(XX+YY-lx-ly+0.1),.0) + np.maximum(100.0*(XX-YY-lx+0.1),.0)                
    V = FunctionSpace(mesh, 'CG', 1)
    phi = Function( V )  
    gdim     = mesh.geometry().dim()
    dofsV    = V.tabulate_dof_coordinates().reshape((-1, gdim))  
    px,py    = [(dofsV[:,0]/lx)*2*Nx, (dofsV[:,1]/ly)*2*Ny]
    dofsV_max =((Nx+1)*(Ny+1) + Nx*Ny)    
    phi = _comp_lsf(px,py,phi,phi_mat,dofsV_max)     

    # Define Omega = {phi<0}     
    class Omega(SubDomain):
        def inside(self, x, on_boundary):
            return .0 <= x[0] <= lx and .0 <= x[1] <= ly and phi(x) < 0

    domains = MeshFunction("size_t", mesh, mesh.topology().dim())
    omega = Omega()
    domains.set_all(0)
    omega.mark(domains, 1)
    dx = Measure('dx')(subdomain_data = domains) 
  
    # Solve elasticity using LU solver
    eps_er, E, nu = [0.001, 1.0, 0.3] 
    mu,lmbda = Constant(E/(2*(1 + nu))),Constant(E*nu/((1+nu)*(1-2*nu)))
    Load = [Point(lx, 0.5)]
    u,v = [TrialFunction(Vvec), TestFunction(Vvec)]
    S1 = 2.0*mu*inner(sym(grad(u)),sym(grad(v))) + lmbda*div(u)*div(v)
    A = assemble( S1*eps_er*dx(0) + S1*dx(1) )   
    b = assemble( inner(Expression(('0.0', '0.0'),degree=2) ,v) * ds(2))    
    U = Function(Vvec)
    delta = PointSource(Vvec.sub(1), Load[0], -1.0)
    delta.apply(b) 
    for bc in bcd: bc.apply(A,b)    

    time1 = time.time()    
    solve(A,U.vector(), b,'lu')  
    time2 =  time.time()
    print('time for solving elasticity system: ', time2-time1)
    
#-----------------------------------------------------------------------
def _comp_lsf(px,py,phi,phi_mat,dofsV_max):               
    for dof in range(0,dofsV_max):              
        if np.rint(px[dof]) %2 == .0: 
            cx,cy = np.int_(np.rint([px[dof]/2,py[dof]/2]))                                            
            phi.vector()[dof] = phi_mat[cy,cx]
        else:
            cx,cy = np.int_(np.floor([px[dof]/2,py[dof]/2]))                      
            phi.vector()[dof] = 0.25*(phi_mat[cy,cx] + phi_mat[cy+1,cx]\
              + phi_mat[cy,cx+1] + phi_mat[cy+1,cx+1])    
    return phi       
    
# ----------------------------------------------------------------------
if __name__ == '__main__':
    _main()
2

Try using MUMPS as the solver rather than UMFPACK.

3

Many thanks. I tried with MUMPS and now I can solve the PDE with large meshes, it solved the “out of memory” problem I had with UMFPACK.

On the other hand using MUMPS the computational time is still around three times larger than when I was using Fenics 2017.

4

If you care about the runtime of your problem, you should consider switching to an iterative solver, as
explained in: How to choose the optimal solver for a PDE problem? - #2 by nate

For elasiticity problems, it is common to use the AMG, see for instance:Bitbucket

5

I’m still baffled why FEniCS 2017 would be “faster” than FEniCS 2019 as measured by the computational linear algebra solve time. How have you installed FEniCS 2019? And most important, how was PETSc installed? Are you sure there’s no other component that’s causing the slowdown?

6

Running your with some slight modificiations, I get the following behavior:

import time
from dolfin import *
import numpy as np
from dolfin import __version__ as version
# ----------------------------------------------------------------------


def _main():

    # Define mesh
    Nx, Ny = [302, 151]   # number of divisions in the x- and y-directions
    lx, ly = [2.0, 1.0]
    XX, YY = np.meshgrid(np.linspace(0.0, lx, Nx+1),
                         np.linspace(0.0, ly, Ny+1))
    # each square of the grid is divided in four triangles
    mesh = RectangleMesh(Point(0.0, 0.0), Point(lx, ly), Nx, Ny, 'crossed')
    Vvec = VectorFunctionSpace(mesh, 'CG', 1)

    # Define boundary conditions
    class DirBd(SubDomain):
        def inside(self, x, on_boundary):
            return near(x[0], .0)
    dirBd = DirBd()
    boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
    boundaries.set_all(0)
    dirBd.mark(boundaries, 1)
    ds = Measure("ds")(subdomain_data=boundaries)
    bcd = [DirichletBC(Vvec, (0.0, 0.0), boundaries, 1)]

    # Define level set function phi
    phi_mat = -np.cos(8.0/lx*pi*XX) * np.cos(4.0*pi*YY) - 0.4\
        + np.maximum(200.0*(0.01-XX**2-(YY-ly/2)**2), .0)\
        + np.maximum(100.0*(XX+YY-lx-ly+0.1), .0) + \
        np.maximum(100.0*(XX-YY-lx+0.1), .0)
    V = FunctionSpace(mesh, 'CG', 1)
    phi = Function(V)
    gdim = mesh.geometry().dim()
    dofsV = V.tabulate_dof_coordinates().reshape((-1, gdim))
    px, py = [(dofsV[:, 0]/lx)*2*Nx, (dofsV[:, 1]/ly)*2*Ny]
    dofsV_max = ((Nx+1)*(Ny+1) + Nx*Ny)
    phi = _comp_lsf(px, py, phi, phi_mat, dofsV_max)

    # Define Omega = {phi<0}
    class Omega(SubDomain):
        def inside(self, x, on_boundary):
            return .0 <= x[0] <= lx and .0 <= x[1] <= ly and phi(x) < 0

    domains = MeshFunction("size_t", mesh, mesh.topology().dim())
    omega = Omega()
    domains.set_all(0)
    omega.mark(domains, 1)
    dx = Measure('dx')(subdomain_data=domains)

    # Solve elasticity using LU solver
    eps_er, E, nu = [0.001, 1.0, 0.3]
    mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1+nu)*(1-2*nu)))
    Load = [Point(lx, 0.5)]
    u, v = [TrialFunction(Vvec), TestFunction(Vvec)]
    S1 = 2.0*mu*inner(sym(grad(u)), sym(grad(v))) + lmbda*div(u)*div(v)
    A = assemble(S1*eps_er*dx(0) + S1*dx(1))
    b = assemble(inner(Expression(('0.0', '0.0'), degree=2), v) * ds(2))
    U = Function(Vvec)
    delta = PointSource(Vvec.sub(1), Load[0], -1.0)
    delta.apply(b)
    for bc in bcd:
        bc.apply(A, b)

    time1 = time.perf_counter()
    solve(A, U.vector(), b, 'lu')
    time2 = time.perf_counter()
    print(
        f'Version: {version} time for solving elasticity system {V.dim()}x{V.dim()}: {time2-time1}')

# -----------------------------------------------------------------------


def _comp_lsf(px, py, phi, phi_mat, dofsV_max):
    for dof in range(0, dofsV_max):
        if np.rint(px[dof]) % 2 == .0:
            cx, cy = np.int_(np.rint([px[dof]/2, py[dof]/2]))
            phi.vector()[dof] = phi_mat[cy, cx]
        else:
            cx, cy = np.int_(np.floor([px[dof]/2, py[dof]/2]))
            phi.vector()[dof] = 0.25*(phi_mat[cy, cx] + phi_mat[cy+1, cx]
                                      + phi_mat[cy, cx+1] + phi_mat[cy+1, cx+1])
    return phi


# ----------------------------------------------------------------------
if __name__ == '__main__':
    _main()

Giving

Version: 2019.2.0.dev0 time for solving elasticity system 91658x91658: 1.1078817000000072

using docker run -ti -v $(pwd):/home/shared -w /home/shared --rm quay.io/fenicsproject/dev:latest
and for 2017.1.0 (docker run -ti -v $(pwd):/home/shared -w /home/shared --rm quay.io/fenicsproject/stable:2017.1.0):
(slight modificiation to print: print("Version: ", version, "time for solving elasticity system ", V.dim(), "x", V.dim(), ": ", time2-time1))

Version:  2017.1.0 time for solving elasticity system  91658 x 91658 :  1.599782442999981
1 Like
7

In 2017 I was running Fenics using Docker. Now I am running Fenics using Conda. Could that be because of Conda?

I just ran the same code on another computer, using UMFPACK and also running Fenics with Conda, and the computational time is similar (around 22s for a 202x101 mesh).

8

Thanks! So you ran these tests using Docker. In 2017 I was also using Docker, now I use Conda to run Fenics. So maybe using Conda is the issue? I will try using Docker.