Iterative solver performance severely reduced using MixedElement

1

In anticipation of adding more equations I’m trying to make use of the MixedElement by default. For example with the ft06 linear elasticity tutorial I simply replace “V = VectorFunctionSpace(mesh, ‘P’, 1)” with

V1 = VectorElement('CG', mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, MixedElement([V1]))

this works, but when the mesh is much denser I notice that iterative solver performance is a magnitude worse (slower) than not using MixedElement, for example using

solve(a == L, u, bc, solver_parameters={"linear_solver": "cg", "preconditioner": "amg"})

It doesn’t seem that this should be the case (as the system to solve should be identical), but I’m possibly missing something fundamental here. Is there some way to still use MixedElement without potentially limiting iterative solver performance?

2

I would suggest moving to dolfinx, as the impact on solver performance on MixedElements have been reduced, as shown in the following minimal example:

import dolfinx
import dolfinx.io
import ufl
from mpi4py import MPI
import numpy as np
import time

def problem(Nx, Ny, Nz, vector_space:bool=True):
    # Scaled variables
    L = 1; W = 0.2
    mu = 1
    rho = 1
    delta = W/L
    gamma = 0.4*delta**2
    beta = 1.25
    lambda_ = beta
    g = gamma

    # Create mesh and define function space
    mesh = dolfinx.BoxMesh(MPI.COMM_WORLD, np.array([[0, 0, 0],[L, W, W]]), [Nx, Ny, Nz], 
    dolfinx.mesh.CellType.tetrahedron)
    if vector_space:
        V = dolfinx.VectorFunctionSpace(mesh, ('P', 1))
    else:
        V1 = ufl.VectorElement('CG', mesh.ufl_cell(), 1)
        V = dolfinx.FunctionSpace(mesh, ufl.MixedElement([V1]))
    # Define boundary condition


    def clamped_boundary(x):
        return np.isclose(x[0], 0)
    if V.num_sub_spaces() == 1:
        V0 = V.sub(0).collapse()
        u_bc = dolfinx.Function(V0)
        u_bc.x.array[:] = 0
        clamped_dofs = dolfinx.fem.locate_dofs_geometrical((V.sub(0),V0), clamped_boundary)
        bc = dolfinx.DirichletBC(u_bc, clamped_dofs, V.sub(0))
    else:
        clamped_dofs = dolfinx.fem.locate_dofs_geometrical(V, clamped_boundary)
        u_bc = dolfinx.Function(V)
        u_bc.x.array[:] = 0
        bc = dolfinx.DirichletBC(u_bc, clamped_dofs)

    # Define strain and stress

    def epsilon(u):
        return 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
        #return sym(nabla_grad(u))

    def sigma(u):
        return lambda_*ufl.nabla_div(u)*ufl.Identity(d) + 2*mu*epsilon(u)

    # Define variational problem
    u = ufl.TrialFunction(V)
    d = len(u)  # space dimension
    v = ufl.TestFunction(V)
    f = dolfinx.Constant(mesh, (0, 0, -rho*g))
    T = dolfinx.Constant(mesh, (0, 0, 0))
    a = ufl.inner(sigma(u), epsilon(v))*ufl.dx
    L = ufl.dot(f, v)*ufl.dx + ufl.dot(T, v)*ufl.ds

    problem = dolfinx.fem.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "cg", "pc_type": "gamg"}) 
    start =time.time()    
    u = problem.solve()
    end = time.time()
    print(f"Vector: {vector_space}, {Nx}x{Ny}x{Nz}: Solve time: {end-start}")
    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, f"u_{Nx}_{Ny}_{Nz}_{vector_space}.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        if V.num_sub_spaces()==1:
            xdmf.write_function(u.sub(0))          
        else:
            xdmf.write_function(u)
    


if __name__ == "__main__":
    Nxs = [10,20,40]
    Nys = [3,6,12]
    Nzs = [3,6,12]
    for Nx, Ny, Nz in zip(Nxs, Nys, Nzs):
       for vector in [True, False]:
           problem(Nx, Ny, Nz, vector)

yielding

Vector: True, 10x3x3: Solve time: 0.0059337615966796875
Vector: False, 10x3x3: Solve time: 0.008495569229125977
Vector: True, 20x6x6: Solve time: 0.05466961860656738
Vector: False, 20x6x6: Solve time: 0.06971216201782227
Vector: True, 40x12x12: Solve time: 0.5019056797027588
Vector: False, 40x12x12: Solve time: 0.5771362781524658

For more information about DOLFINx, you could have a look at my tutorial: The FEniCSx tutorial — FEniCSx tutorial

One should also distinguish between solver performance and the time it takes to assemble the left and right hand sides, and actually run the iterative solver.

1 Like