Crack growth rate validation

This example demonstrates how to compute the rate of advancing the crack using the asymptotic approach, and compares with Monte Carlo calculations.

Isometric ensemble

This is first done in the isometric thermodynamic ensemble. The crack growth rate will be computed for a set of increasing nondimensional bond energies:

varepsilons = [10, 25, 100, 1000]

Now the rate (relative to the initial rate) is computed as a function of the nondimensional end separation in a loop for each model. The default calculation approach is asymptotic:

import numpy as np
from statmechcrack import Crack

v = np.linspace(1, 11, 33)
k_asymptotic = np.zeros((len(v), len(varepsilons)))
for i, varepsilon in enumerate(varepsilons):
    model = Crack(varepsilon=varepsilon)
    k_asymptotic[:, i] = model.k(v, ensemble='isometric')

Monte Carlo is another calculation approach and does not make approximations, but takes a considerable amount of time. These calculations have already been performed, and the results archived in the local data folder. The steps to obtain these results are commented below, and the results are simply loaded instead:

k_monte_carlo = np.zeros((len(v), len(varepsilons)))
for i, varepsilon in enumerate(varepsilons):
#     model = Crack(varepsilon=varepsilon)
#     k_monte_carlo[:, i] = model.k(
#         v, ensemble='isometric', approach='monte carlo',
#         num_processes=32, num_samples=int(3.125e7))
    k_monte_carlo[:, i] = np.load(
        'data/k_monte_carlo_isometric_varepsilon_'
        + str(varepsilon) + '.npy')[2]

The results can now be plotted together:

import matplotlib.pyplot as plt

_ = plt.semilogy(v, k_asymptotic)
_ = plt.semilogy(v, k_monte_carlo, 'k--')

The relative L2 error can be plotted across more nondimensional energies:

varepsilons = [10, 15, 25, 40, 65, 100, 160, 250, 400, 630, 1000, 1600, 2500]
relative_error = np.zeros(len(varepsilons))
for i, varepsilon in enumerate(varepsilons):
    k_monte_carlo = np.load(
        'data/k_monte_carlo_isometric_varepsilon_'
        + str(varepsilon) + '.npy')[2]
    v = np.load(
        'data/k_monte_carlo_isometric_varepsilon_'
        + str(varepsilon) + '.npy')[0]
    k_asymptotic = Crack(varepsilon=varepsilon).k(v, ensemble='isometric')
    _ = plt.semilogy(v, k_asymptotic)
    _ = plt.semilogy(v, k_monte_carlo, 'k--')
    relative_error[i] = np.sqrt(
        np.trapz((np.log(k_asymptotic/k_monte_carlo))**2, x=v) /
        np.trapz(np.log(k_monte_carlo)**2, x=v))
plt.show()
_ = plt.loglog(varepsilons, relative_error, 'o-')

Isotensional ensemble

varepsilons = [10, 25, 100, 1000]
v = np.linspace(1, 11, 33)
p = np.zeros((len(v), len(varepsilons)))
k_asymptotic = np.zeros((len(v), len(varepsilons)))
k_monte_carlo = np.zeros((len(v), len(varepsilons)))
for i, varepsilon in enumerate(varepsilons):
    model = Crack(varepsilon=varepsilon)
    p[:, i] = model.p(v, ensemble='isometric')
    k_asymptotic[:, i] = model.k(p[:, i], ensemble='isotensional')
#     model = Crack(varepsilon=varepsilon)
#     k_monte_carlo[:, i] = model.k(
#         v, ensemble='isotensional', approach='monte carlo',
#         num_processes=32, num_samples=int(3.125e7))
    k_monte_carlo[:, i] = np.load(
        'data/k_monte_carlo_isotensional_varepsilon_'
        + str(varepsilon) + '.npy')[2]
_ = plt.semilogy(p, k_asymptotic)
_ = plt.semilogy(p, k_monte_carlo, 'k--')

varepsilons = [10, 15, 25, 40, 65, 100, 160, 250, 400, 630, 1000, 1600, 2500]
relative_error = np.zeros(len(varepsilons))
for i, varepsilon in enumerate(varepsilons):
    k_monte_carlo = np.load(
        'data/k_monte_carlo_isotensional_varepsilon_'
        + str(varepsilon) + '.npy')[2]
    p = np.load(
        'data/k_monte_carlo_isotensional_varepsilon_'
        + str(varepsilon) + '.npy')[0]
    k_asymptotic = Crack(varepsilon=varepsilon).k(p, ensemble='isotensional')
    _ = plt.semilogy(p, k_asymptotic)
    _ = plt.semilogy(p, k_monte_carlo, 'k--')
    relative_error[i] = np.sqrt(
        np.trapz((np.log(k_asymptotic/k_monte_carlo))**2, x=v) /
        np.trapz(np.log(k_monte_carlo)**2, x=v))
plt.show()
_ = plt.loglog(varepsilons, relative_error, 'o-')