Python tutorial - GeoEnergyLab-EPFL/BigWham GitHub Wiki

Penny shaped crack (Axi symmteric kernel)

from bigwham4py import BEMatrix
import numpy as np

radius = 1.0
pressure = 1.0

G = 1.0
nu = 0.0
E = (2 * G) * (1 + nu)

max_leaf_size = 20
eta = 5.0
eps_aca = 1e-4

nelts = 1000
coor1D = np.linspace(0, radius, nelts + 1)
coor = np.transpose(np.array([coor1D, coor1D * 0.0]))
conn = np.fromfunction(lambda i, j: i + j, (nelts, 2), dtype=np.int_)

# Create BE-Matrix
kernel = "Axi3DS0-H"
hmat = Hmatrix(kernel, coor, conn, np.array([E, nu]), max_leaf_size, eta, eps_aca)

col_pts = hmat.getMeshCollocationPoints()

pre_fac = (8 * (1 - nu * nu)) / (np.pi * E)
dd = np.zeros(col_pts.shape)
dd[:, 1] = pre_fac * np.sqrt(radius * radius - np.linalg.norm(col_pts[:, :], axis=1)**2)


# calculate tractions
t = hmat.matvec(dd.flatten())

# print(dd)
# print(t)

t_anal = np.zeros(col_pts.shape)
t_anal[:, 1] = pressure

print("L2 Rel error {}".format(np.linalg.norm(t - t_anal.flatten()) / nelts))


# %% Plot DD
r = np.linalg.norm(col_pts[:, :], axis=1)
plt.figure()
plt.plot(r, dd[:, :], ".")
plt.xlabel("radius")
plt.ylabel("DD")


# %% Plot traction
plt.figure()
plt.plot(r, t.reshape(-1, 2)[:, :], ".")
plt.xlabel("radius")
plt.ylabel("traction")


# %% Plot traction error
plt.figure()
plt.plot(r, np.abs(t.reshape(-1, 2)[:, :] - t_anal) / pressure, ".")
plt.semilogy()
plt.xlabel("radius")
plt.ylabel("traction error")

Scaling with OMP threads for Penny shaped crack (Axi symmteric kernel)

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# %%
import os
import sys
import matplotlib.pyplot as plt

print(sys.path)

#%%


home = os.environ["HOME"]

sys.path.append(home + "/geolab/BigWham/build/interfaces/python")


from bigwham4py import BEMatrix
import numpy as np





radius = 1.0
pressure = 1.0

G = 1.0
nu = 0.0
E = (2 * G) * (1 + nu)


# %% Mesh Information
nelts = 2000
coor1D = np.linspace(0, radius, nelts + 1)
coor = np.transpose(np.array([coor1D, coor1D * 0.0]))
conn = np.fromfunction(lambda i, j: i + j, (nelts, 2), dtype=np.int_)

# %% Create H-matrix
max_leaf_size = 100
eta = 3.0
eps_aca = 1e-4
kernel = "Axi3DS0-H"

hmat = BEMatrix(kernel, coor, conn, np.array([E, nu]), max_leaf_size, eta, eps_aca)

# %% Compression ratio
cr = hmat.getCompression()
print("Compression ratio: ", cr)
# print("Number of NNZs: ", hsparse.nnz)

# %% Calculate DD
col_pts = hmat.getMeshCollocationPoints()
pre_fac = (8 * (1 - nu * nu)) / (np.pi * E)
dd = np.zeros(col_pts.shape)
dd[:, 0] = pre_fac * np.sqrt(
    radius * radius - np.linalg.norm(col_pts[:, :], axis=1) ** 2
)
dd_fl = dd.flatten()


# %% Measure time with varying number of threads
from threadpoolctl import threadpool_limits
import time
ntimes = 10000
for threads in [2, 4, 8, 16, 32, 64, 96, 128]:  # Varying number of threads
    with threadpool_limits(limits=threads, user_api='openmp'):
        hmat.get_omp_threads()
        hmat = BEMatrix(kernel, coor, conn, np.array([E, nu]), max_leaf_size, eta, eps_aca)
        t1 = time.time()
        for i in range(ntimes):
            hmat @ dd_fl
        elapsed = (time.time() - t1)
        print(f"Threads: {threads}, Elapsed Time: {elapsed:.4e} seconds")