<jupyter_start><jupyter_text>Measuring Properties of Ball Spring SolidTh<NAME><jupyter_code>from IPython.display import HTML

def playVideo(path):
    return HTML("""
                <video width="320" height="240" controls>
                  <source src="{}" type="video/mp4">
                </video>
                """.format(path))<jupyter_output><empty_output><jupyter_text>What do we want to measure Youngs Modulus    - A measure of the stiffness of a material $$E = \frac{Stress}{Strain}$$ $$Stress = \frac{A}{F} \;\;\;\;\;\;\;\;\;\;\;\;\;\;Strain = \frac{\Delta L}{L}$$ Coefficient of Restitution    - A measure of kinetic energy lost in collision $$e = \sqrt{\frac{T_{i}}{T_{f}}}$$<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mplEasyAnimate import animation
from tqdm import tqdm
from scipy.special import comb
from scipy.special import factorial
from numba import njit, jit
@njit
def rk4(y0, IDx, IDy, h, sep, pList, massList, ts, b, k, F):
    k1 = h*nbody(y0, IDx, IDy, h, sep, pList, massList, ts, b, k, F)
    k2 = h*nbody(y0+k1/2, IDx, IDy, h, sep, pList, massList, ts, b, k, F)
    k3 = h*nbody(y0+k2/2, IDx, IDy, h, sep, pList, massList, ts, b, k, F)
    k4 = h*nbody(y0+k3, IDx, IDy, h, sep, pList, massList, ts, b, k, F)
    return y0 + (k1/6)+(k2/3)+(k3/3)+(k4/6)
@njit
def getVecs(r1, r2):
    sr = r2-r1
    srmag = np.sqrt(sr[0]**2+sr[1]**2)
    return sr, srmag, sr/srmag
@njit
def getForceComponents(r, v, pList, IDx, IDy, k, b, sep, xoffsets, yoffsets, ts):
    FMag = np.zeros(shape=(2,))
    for oID in range(len(xoffsets)):
        dx = xoffsets[oID]
        dy = yoffsets[oID]
        sr, srmag, srhat = getVecs(r, pList[IDx+dx, IDy+dy, :2])
        if srmag > 0.1:
            FMag += -k*(sep-srmag)*srhat-b*v
    return FMag<jupyter_output><empty_output><jupyter_text>Model  Boundary ConditionsThere are eight boundary conditions we need to handel seperatly     1) Corners (4x)    2) Edges (4x)    3) Remainder Corners<jupyter_code>if IDx == IDy == 0:
        FMag = getForceComponents(r, v, pList, IDx, IDy, k, b, sep, [1, 0], [0, 1], ts)
        FMag = np.array([0.0, 0.0])

    elif IDx == pList.shape[0]-1 and IDy == pList.shape[1]-1:
        FMag = getForceComponents(r, v, pList, IDx, IDy, k, b, sep, [-1, 0], [0, -1], ts) + F
        
    elif IDx == 0 and IDy == pList.shape[1]-1:
        FMag = getForceComponents(r, v, pList, IDx, IDy, k, b, sep, [1, 0], [0, -1], ts)
        FMag = np.array([0.0, 0.0])

    elif IDx == pList.shape[0]-1 and IDy == 0:
        FMag = getForceComponents(r, v, pList, IDx, IDy, k, b, sep, [-1, 0], [0, 1], ts) + F<jupyter_output><empty_output><jupyter_text>Edges<jupyter_code>elif IDy == 0 and IDx != 0 and IDx != pList.shape[0]-1:
        FMag = getForceComponents(r, v, pList, IDx, IDy, k, b, sep, [-1, 0, 1], [0, 1, 0], ts)

    elif IDy == pList.shape[1]-1 and IDx != 0 and IDx != pList.shape[0]-1:
        FMag = getForceComponents(r, v, pList, IDx, IDy, k, b, sep, [-1, 0, 1], [0, -1, 0], ts)

    elif IDx == 0 and IDy != 0 and IDy != pList.shape[1]-1:
        FMag = getForceComponents(r, v, pList, IDx, IDy, k, b, sep, [0, 1, 0], [-1, 0, 1], ts)
        FMag = np.array([0.0, 0.0])
        
    elif IDx == pList.shape[0]-1 and IDy != 0 and IDy != pList.shape[1]-1:
        FMag = getForceComponents(r, v, pList, IDx, IDy, k, b, sep, [0, -1, 0], [-1, 0, 1], ts) + F<jupyter_output><empty_output><jupyter_text>All Other Particles<jupyter_code>else:
        FMag = getForceComponents(r, v, pList, IDx, IDy, k, b, sep, [0, 1, 0, -1], [-1, 0, 1, 0], ts)
def int_n_model(pList, mass, sep, h, tf=1, b=0.2, k=10, F=np.array([0,0]), pbar=True):
    ts = np.arange(0, tf, h)
    ys = np.zeros(shape=(len(ts)+1, pList.shape[0], pList.shape[1], pList.shape[2]))
    ys[0] = pList

    for i in tqdm(range(ts.shape[0]), disable=not pbar):
        for IDx in range(ys[i, :].shape[0]):
            for IDy in range(ys[i, :].shape[1]):
                ys[i+1, IDx, IDy] = rk4(ys[i, IDx, IDy], IDx, IDy, h, sep, ys[i, :, :], massList, i, b, k, F)
    return np.arange(0, tf+h, h), ys
def make_lattice(Nv, Nh, sep):
    lattice = np.zeros(shape=(Nv, Nh, 4))
    for nv in range(Nv):
        for nh in range(Nh):
            lattice[nh, nv] = np.array([nh*sep, nv*sep, 0, 0])
    
    return lattice
def plot_system(y, F=None):
    fig = plt.figure(figsize=(10, 7))
    ax = fig.add_subplot(111)
    ax.plot(y[:, 1:-1, 0], y[:, 1:-1, 1], 'C0o')
    if F:
        ax.quiver(y[-1, 1:-1, 0], y[-1, 1:-1, 1], F[0], F[1])
    return fig, ax
def animate_system(filename, ys, F=None, skip=1):
    anim = animation(filename, fps=30)
    for i, y in tqdm(enumerate(ys), total=ys.shape[0]):
        if i%skip == 0:
            fig, ax = plot_system(y, F=F)
            ax.set_xlim(np.min(ys[:, 0, 0])-0.1, np.max(ys[:, -1, 0])+0.1)
            ax.set_ylim(np.min(ys[0, :, 1])-0.1, np.max(ys[-1, :, 1])+0.1)
            anim.add_frame(fig)
            plt.close(fig)
    
    anim.close()<jupyter_output><empty_output><jupyter_text>Initial ContidionsI defined helper funcions to build a latice of particles seperated by the equilibrium length some spring I connecte them with<jupyter_code>sep=1
n = 10
lattice = make_lattice(n, n, sep)
massList = np.ones(shape=(n, n))
F = np.array([0, 0])
fig, ax = plot_system(lattice)
lattice[9, :, 0] += 1
fig, ax = plot_system(lattice)
ts, ys = int_n_model(lattice, massList, sep, 0.01, tf=15, b=1, k=10)
animate_system('PresentationAnimations/StretchAndRelease.mp4', ys, skip=5)
playVideo('PresentationAnimations/StretchAndRelease.mp4')
fig = plt.figure(figsize=(10, 7))
plt.plot(ts, ys[:, -1, 0, 0])
plt.xlabel('Time [s]', fontsize=17)
plt.ylabel('X Pos [m]', fontsize=17)
plt.savefig('PresentationFigures/StretchAndReleaseXpos.png', dpi=350, bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>This is an overdamped system Calculation of Youngs ModulusRecall that: $$E = \frac{Stress}{Strain}$$ $$Stress = \frac{A}{F} \;\;\;\;\;\;\;\;\;\;\;\;\;\;Strain = \frac{\Delta L}{L}$$ MethodWe want to measure youngs modulus for the material and see how it varies with spring constant    1) Pin one end of the solid    2) Apply a constant force to the other over some integration with a given spring constant    3) Measure the change in length and calculate the strain to get Youngs Modulus    4) Repeat 2-3 for a set of spring constants<jupyter_code>playVideo('PresentationAnimations/ConstantForce.mp4')
def get_strain(ys):
    Li = np.mean(ys[0, -1, :, 0]-ys[0, 0, :, 0])
    Lf = np.mean(ys[-1, -1, :, 0]-ys[-1, 0, :, 0])
    dL = Lf-Li
    return dL/Li
def get_stress(ys, F):
    A = ys[-1, -1, -1, 1]-ys[-1, -1, 0, 1]
    FMag = np.sqrt(F[0]**2+F[1]**2)
    return FMag/A
def get_youngs_modulus(ys, F):
    return get_stress(ys, F)/get_strain(ys)
ea = list()
K = np.linspace(1, 20, 100)
F = np.array([10, 0])
for k in tqdm(K):
    ts, ys = int_n_model(lattice, massList, sep, 0.01, tf=15, b=1, k=k, F=F, pbar=False)
    ea.append(get_youngs_modulus(ys, F))
plt.figure(figsize=(10, 7))
plt.plot(np.linspace(1, 20, 100), ea)
plt.xlabel(r'Spring Constant [N m$^{-1}$]', fontsize=17)
plt.ylabel(r'Youngs Modulus [kg m$^{-2}$ s$^{-2}$]', fontsize=17)
plt.savefig('PresentationFigures/YoungsModulusVsK.png', dpi=350, bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>Set of spring constants<jupyter_code>for k in tqdm(np.linspace(1, 20, 100)):
    ts, ys = int_n_model(lattice, massList, sep, 0.01, tf=15, b=1, k=k, F=F, pbar=False)
    e.append(get_youngs_modulus(ys, F))<jupyter_output><empty_output><jupyter_text>Coefficient of RestitutionNo results yetWill reverse the momentum of the bottom row of springs when the entire solid has some velocity Honors Projectnbody simulations are naievly an $O(n^{2})$ calculationspython is slow when running for loops Ways around this        1) Ignore it, just use pure Python    2) Use numpy for a small speed increase    3) Make your code a little more C like and use numba    4) Use pycuda    5) just give in, rewrite it in C, and the pure CUDA, and       then drop out of College, I bet Home Depot is Hiring. I never used pure python, but I will make a speed comparision between 2-5<jupyter_code>times = {'Numpy': 723, 'Numba': 20, 'PyCUDA': 1.6, 'C': 1.52, 'CUDA': 0.59}
fig = plt.figure(figsize=(10, 7))
plt.bar(range(len(times)), list(times.values()), align='center')
plt.xticks(range(len(times)), list(times.keys()))
plt.ylabel('Sys Time to complete 1000 time steps\nfor 100 particles [s]', fontsize=20)
plt.savefig('PresentationFigures/TimesWithNumpy.png', dpi=350, bbox_inches='tight')
times = {'Numba': 20, 'PyCUDA': 1.6, 'C': 1.52, 'CUDA': 0.59}
fig = plt.figure(figsize=(10, 7))
plt.bar(range(len(times)), list(times.values()), align='center')
plt.xticks(range(len(times)), list(times.keys()))
plt.ylabel('Sys Time to complete 1000 time steps\nfor 100 particles [s]', fontsize=20)
plt.savefig('PresentationFigures/TimesWithNumba.png', dpi=350, bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1768617021276596, 'fraction_numerical': 0.041555851063829786, 'mean_word_length': 3.659266907589055, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '383510', 'license': 'MIT', 'n_tokens_mistral': 3919, 'n_tokens_neox': 3661, 'n_words': 1000, 'path': 'YoungsModulusOfABallSpringSolid/Presentation.ipynb', 'repo_name': 'tboudreaux/ComputationalPhysics', 'size': '9035'}