Assembly times in serial vs parallel

[bhaveshshrimali]

Following the discussion in #587 (comment) I tried to do some benchmarking on my laptop Intel(R) Xeon(R) CPU E3-1535M v6 @ 3.10GHz with 4 cores

and got some rather interesting numbers.. (performance is outstanding, but the differences in serial vs parallel are unclear to me)

parallel=False

# Assembly Runtimes

# Assembly Runtimes

|   DOF   |             FElupe               |           scikit-fem            |
| ------- | -------------------------------- | ------------------------------- |
|     24  |                0.01              |                 nan             |
|    192  |                0.01              |                 nan             |
|   1536  |                0.05              |                 nan             |
|  12288  |                0.30              |                 nan             |
|  98304  |                1.95              |                 nan             |

| Analysis       |        DOF/s      |
| -------------- | ----------------- |
| FElupe         |  30330 +/-15814   |
| skfem          |     nan +/-  nan  |

parallel = True

# Assembly Runtimes

|   DOF   |             FElupe               |           scikit-fem            |
| ------- | -------------------------------- | ------------------------------- |
|     24  |                0.01              |                 nan             |
|    192  |                0.02              |                 nan             |
|   1536  |                0.04              |                 nan             |
|  12288  |                0.27              |                 nan             |
|  98304  |                2.18              |                 nan             |

| Analysis       |        DOF/s      |
| -------------- | ----------------- |
| FElupe         |  28130 +/-17791   |
| skfem          |     nan +/-  nan  |

Code (you can ignore non-standard imports for scikit-fem)

# %%
from timeit import timeit

import matplotlib.pyplot as plt
import numpy as np
from meshio import xdmf
import felupe as fem
from einsumt import einsumt
from felupe.math import det, inv, transpose, identity

# from skfem import (
#     MeshTet,
#     ElementVectorH1,
#     ElementTetP2,
#     ElementVector,
#     ElementTetP1,
#     Basis,
#     asm,
# )
from scipy.sparse import bmat
from warnings import filterwarnings

import sys

sys.path.append("../")
# from scikitfemEx36Numba import F1, F2, A11, A22, A12, a1, a2, b11, b12, b22
# from skfem.io import from_meshio

filterwarnings("ignore")

parallel = True

import logging

logging.basicConfig(
    filename=f"output_par_{str(parallel)}.log",
    level=logging.INFO,
    format="%(message)s",
)


def log(message):
    logging.info(message)


# %%
def load_mesh(n):
    # mesh = MeshTet().refined(n)
    # mesh.save(f"mesh_{n:d}.xdmf")
    # meshio_mesh = xdmf.read(f"mesh_{n:d}.xdmf")
    # mesh_felupe = fem.Mesh(
    #     points=meshio_mesh.points,
    #     cells=meshio_mesh.get_cells_type("tetra"),
    #     cell_type=meshio_mesh.cells[0].type,
    # )
    # region = fem.RegionTetra(mesh=mesh_felupe)
    # mesh_felupe = mesh_felupe.flip(np.any(region.dV < 0, axis=0))
    mesh_felupe = fem.Cube(n=n).triangulate()
    mesh_felupe.save(f"mesh_{n:d}.xdmf")
    return mesh_felupe


def pre_hyperelastic(n, **kwargs):
    mesh = load_mesh(n)
#    mesh = mesh.add_midpoints_edges()
    region = fem.RegionTetra(mesh)
    field = fem.FieldContainer([fem.Field(region, dim=3)])
    umat = fem.NeoHooke(mu=1.0, bulk=1.0e4)
    solid = fem.SolidBody(umat, field)
    return mesh, solid


def pre_skfem_hyperelastic(n):
    mesh = from_meshio(xdmf.read(f"mesh_{n:d}.xdmf"))

    uelem = ElementVector(ElementTetP2())
    pelem = ElementTetP1()
    elems = {"u": uelem, "p": pelem}
    basis = {field: Basis(mesh, e, intorder=1) for field, e in elems.items()}

    du = basis["u"].zeros()
    dp = basis["p"].zeros()
    uv = basis["u"].interpolate(du)
    pv = basis["p"].interpolate(dp)
    # t1 = time()
    # K11 = asm(b11, basis["u"], basis["u"], disp=uv, press=pv)
    # K12 = asm(b12, basis["p"], basis["u"], disp=uv, press=pv)
    # K22 = asm(b22, basis["p"], basis["p"], disp=uv, press=pv)
    # K = bmat([[K11, K12], [K12.T, K22]], "csr")
    return (b11, b12, b22), basis, uv, pv


def assemble_matrix_skfem(forms, basis, uv, pv):
    b11, b12, b22 = forms
    K11 = asm(b11, basis["u"], basis["u"], disp=uv, press=pv)
    K12 = asm(b12, basis["p"], basis["u"], disp=uv, press=pv)
    K22 = asm(b22, basis["p"], basis["p"], disp=uv, press=pv)
    K = bmat([[K11, K12], [K12.T, K22]], "csr")
    return K


# %%
# pre_hyperelastic(3)

# %%

log("# Assembly Runtimes")
log("")
log("|   DOF   |             FElupe               |           scikit-fem            |")
log("| ------- | -------------------------------- | ------------------------------- |")

points_per_axis = np.array([2, 4, 8, 16, 32])

number = 3


runtimes = np.zeros((len(points_per_axis), 2))

for i, n in enumerate(points_per_axis):
    mesh, solid = pre_hyperelastic(n)
    matrix = solid.assemble.matrix(parallel=parallel)
    time_hyperelastic = (
        timeit(lambda: solid.assemble.matrix(parallel=parallel), number=number) / number
    )
    # forms, basis, uv, pv = pre_skfem_hyperelastic(n)
    # time_hyperelastic_skfem = (
    #     timeit(lambda: assemble_matrix_skfem(forms, basis, uv, pv), number=number)
    #     / number
    # )
    runtimes[i] = time_hyperelastic, np.nan

    log(
        f"|{mesh.points.size:7d}  | {runtimes[i][0]:19.2f}              | {runtimes[i][1]:19.2f}             |"
    )

dofs_le = points_per_axis**3 * 3 / runtimes[:, 0]
dofs_he = points_per_axis**3 * 3 / runtimes[:, 1]

log("")
log("| Analysis       |        DOF/s      |")
log("| -------------- | ----------------- |")
log(f"| FElupe         |  {np.mean(dofs_le):5.0f} +/-{np.std(dofs_le):5.0f}   |")
log(f"| skfem          |   {np.mean(dofs_he):5.0f} +/-{np.std(dofs_he):5.0f}  |")


# %%
plt.figure()
# plt.loglog(
#     points_per_axis**3 * 3,
#     runtimes[:, 1],
#     "C0",
#     label=r"skfem",
# )
plt.loglog(
    points_per_axis**3 * 3,
    runtimes[:, 0],
    "C1--",
    label=r"FElupe",
)
plt.xlabel(r"Number of degrees of freedom $\longrightarrow$")
plt.ylabel(r"Runtime in s $\longrightarrow$")
plt.legend()
plt.tight_layout()
# plt.savefig("benchmark.svg", transparent=False)
plt.savefig(f"benchmark_par_{str(parallel)}.png", transparent=False)

[adtzlr]

Hi @bhaveshshrimali!

Yes, thanks for pointing that out! I got similar results on my Windows machine for RegionTetra. However, for RegionHexahedron the parallel keyword improves performance already by about 1.5 and for RegionQuadraticHexahedron by about 2.5. For small meshes, parallel introduces some extra overhead.

[adtzlr]

In older versions of FElupe the difference was bigger. Recently, I tried to enhance the default parallel=False performance by using the out-keyword in numpy-functions when possible. Hopefully this was successful 😃

[adtzlr]

I'm not sure why this is happening actually. In RegionTetra, there is only one quadrature point. Nothing to sum here on integration compared to Hexahedrons. That is the only difference from a data point of view.

[bhaveshshrimali]

Yeah, even I was thrown away with the numbers initially. I tried to compare both Tet and Hex assemblies and noticed similar times (though not strictly benefitting in parallel)

parallel=False

DOF (hex)	DOF (tet)	FElupe (hex)	FElupe (tet)
24	24	0.00572423	0.00548673
192	192	0.0103467	0.0079608
1536	1536	0.0516957	0.0347423
12288	12288	0.375947	0.252385
98304	98304	3.78216	2.00618
786432	786432	33.3451	17.2463

parallel=True

DOF (hex)	DOF (tet)	FElupe (hex)	FElupe (tet)
24	24	0.0079882	0.00683567
192	192	0.0126934	0.0106748
1536	1536	0.0472516	0.0401419
12288	12288	0.31048	0.248962
98304	98304	3.19302	2.45148
786432	786432	29.4799	19.3428

While, hex improves when going from setting parallel=False to True, not so much for Tet.

[adtzlr]

#680 could also be interesting - but it is far away from ready to be merged.