{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":39763,"databundleVersionId":11756775,"sourceType":"competition"}],"dockerImageVersionId":30918,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-07T07:58:37.554858Z","iopub.execute_input":"2025-06-07T07:58:37.555308Z","iopub.status.idle":"2025-06-07T07:58:39.065325Z","shell.execute_reply.started":"2025-06-07T07:58:37.555274Z","shell.execute_reply":"2025-06-07T07:58:39.064078Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"This notebook presents a modified version of the seismic wave simulation code from [@jaewook704's notebook](https://www.kaggle.com/code/jaewook704/waveform-inversion-vel-to-seis).\n\nMy goal was to get the Python simulation results to more closely match the original MATLAB code. After these modifications, the numerical difference is now down to about $1 \\times 10^{-5}$ across different velocity models.\n\nFull disclosure: These modifications were developed with the help of Google's Gemini. I honestly just asked it to \"reduce the numerical calculation error between Python and Matlab,\" and while I'm still exploring the theory behind it, the results work very well.","metadata":{}},{"cell_type":"code","source":"# https://arxiv.org/pdf/2111.02926\n# https://csim.kaust.edu.sa/files/SeismicInversion/Chapter.FD/lab.FD2.8/lab.html\n\n# no modification\ndef ricker(f, dt, nt=None):\n    nw = int(2.2 / f / dt)\n    nw = 2 * (nw // 2) + 1\n    nc = nw // 2 + 1  # 중심 인덱스를 1-based 기준으로 설정\n\n    k = np.arange(1, nw + 1)  # 1-based index\n    alpha = (nc - k) * f * dt * np.pi\n    beta = alpha ** 2\n    w0 = (1.0 - 2.0 * beta) * np.exp(-beta)\n\n    # 1-based wavelet 생성\n    if nt is not None:\n        if nt < len(w0):\n            raise ValueError(\"nt is smaller than condition!\")\n        w = np.zeros(nt + 1)  # dummy 포함\n        w[1:len(w0) + 1] = w0\n    else:\n        w = np.zeros(len(w0) + 1)\n        w[1:] = w0\n\n    # 1-based time axis 생성\n    if nt is not None:\n        tw = np.arange(1, len(w)) * dt\n    else:\n        tw = np.arange(1, len(w)) * dt\n\n    return w, tw\n\ndef AbcCoef2D(vel, nbc, dx):\n    \"\"\"\n    Calculates coefficients for a 2D Absorbing Boundary Condition (ABC).\n    This is a Python/NumPy translation of the provided MATLAB function.\n\n    Args:\n        vel (np.ndarray): The padded 2D velocity model.\n        nbc (int): The number of padding cells (boundary width).\n        dx (float): The spatial grid interval.\n\n    Returns:\n        np.ndarray: A 2D array of damping coefficients.\n    \"\"\"\n    nzbc, nxbc = vel.shape\n    velmin = np.min(vel)\n    nz = nzbc - 2 * nbc\n    nx = nxbc - 2 * nbc\n\n    if nbc <= 1:\n        return np.zeros_like(vel)\n\n    a = (nbc - 1) * dx\n    kappa = 3.0 * velmin * np.log(1e7) / (2.0 * a)\n\n    damp1d = kappa * ((np.arange(nbc) * dx / a) ** 2)\n    damp = np.zeros((nzbc, nxbc), dtype=np.float64)\n\n    # Fill left and right damping zones\n    for iz in range(nzbc):\n        damp[iz, 0:nbc] = damp1d[::-1]\n        damp[iz, nx + nbc : nx + 2 * nbc] = damp1d\n\n    # Fill top and bottom damping zones\n    for ix in range(nbc, nbc + nx):\n        damp[0:nbc, ix] = damp1d[::-1]\n        damp[nbc + nz : nz + 2 * nbc, ix] = damp1d\n        \n    return damp\n\ndef padvel(v0, nbc):\n    \"\"\"\n    Pads the velocity model by extending the edge values outward.\n    \"\"\"\n    return np.pad(v0, pad_width=nbc, mode='edge')\n\ndef expand_source(s0, nt):\n    \"\"\"\n    Ensures the source time function has length 'nt'.\n    \"\"\"\n    nt0 = s0.size\n    if nt0 < nt:\n        s = np.zeros(nt, dtype=np.float64)\n        s[:nt0] = s0\n        return s\n    else:\n        return s0[:nt].astype(np.float64)\n\ndef adjust_sr(coord, dx, nbc):\n    \"\"\"\n    Converts physical source/receiver coordinates to grid indices.\n    \"\"\"\n    # MATLAB's round(x.5) rounds away from zero. NumPy's np.round(x.5)\n    # rounds to the nearest even integer. Using np.floor(x + 0.5) for\n    # positive numbers emulates MATLAB's behavior.\n    round_to_int = lambda x: np.floor(x + 0.5).astype(int)\n\n    isx = round_to_int(coord['sx'] / dx) + nbc\n    isz = round_to_int(coord['sz'] / dx) + nbc\n    igx = round_to_int(coord['gx'] / dx) + nbc\n    igz = round_to_int(coord['gz'] / dx) + nbc\n\n    if np.abs(coord['sz']) < 0.5:\n        isz += 1\n        \n    igz += (np.abs(coord['gz']) < 0.5).astype(int)\n    \n    return isx, isz, igx, igz\n\ndef a2d_mod_abc24(v, nbc, dx, nt, dt, s, coord, isFS):\n    \"\"\"\n    Performs a 2D acoustic wave finite-difference simulation (4th order).\n    \"\"\"\n    n_receivers = coord['gx'].size\n    seis = np.zeros((nt, n_receivers), dtype=np.float64)\n    \n    c1 = -2.5\n    c2 = 4.0 / 3.0\n    c3 = -1.0 / 12.0\n    \n    v_padded = padvel(v, nbc)\n    abc = AbcCoef2D(v_padded, nbc, dx)\n    \n    alpha = (v_padded * dt / dx)**2\n    kappa = abc * dt\n    temp1 = 2 + 2 * c1 * alpha - kappa\n    temp2 = 1 - kappa\n    beta_dt = (v_padded * dt)**2\n    \n    s = expand_source(s, nt)\n    isx, isz, igx, igz = adjust_sr(coord, dx, nbc)\n\n    p0 = np.zeros_like(v_padded, dtype=np.float64)\n    p1 = np.zeros_like(v_padded, dtype=np.float64)\n    \n    for it in range(nt):\n        laplacian = (\n            c2 * (np.roll(p1, 1, axis=1) + np.roll(p1, -1, axis=1) +\n                  np.roll(p1, 1, axis=0) + np.roll(p1, -1, axis=0)) +\n            c3 * (np.roll(p1, 2, axis=1) + np.roll(p1, -2, axis=1) +\n                  np.roll(p1, 2, axis=0) + np.roll(p1, -2, axis=0))\n        )\n        \n        p = temp1 * p1 - temp2 * p0 + alpha * laplacian\n        p[isz, isx] += beta_dt[isz, isx] * s[it]\n        \n        if isFS:\n            p[nbc, :] = 0.0\n            p[nbc-1, :] = -p[nbc+1, :]\n            p[nbc-2, :] = -p[nbc+2, :]\n\n        seis[it, :] = p[igz, igx]\n        \n        p0 = p1.copy()\n        p1 = p.copy()\n        \n    return seis\n    \ndef vel_to_seis(vel, method='abc24'):\n    \"\"\"\n    Runs the simulation for multiple sources and collects the seismograms.\n    \n    Args:\n        vel (np.ndarray): The (70, 70) velocity model.\n        method (str): The simulation function to use ('abc24').\n        \n    Returns:\n        np.ndarray: Stacked seismogram data of shape (5, 1001, 70).\n    \"\"\"\n    # 1. Model and Simulation Parameters\n    nz, nx = vel.shape\n    dx = 10.0\n    nbc = 120\n    nt = 1001\n    dt = 1e-3\n    freq = 15.0\n    isFS = False  # Use free surface condition or not\n\n    # 2. Generate Ricker wavelet source\n    s, _ = ricker(freq, dt)\n    \n    # 3. Setup Receiver Coordinates\n    # Receivers are placed at every grid point horizontally at a fixed depth.\n    coord = {}\n    coord['sz'] = 1 * dx\n    coord['gx'] = np.arange(nx) * dx\n    coord['gz'] = np.ones(nx) * dx\n    \n    # 4. Loop over source positions and run simulation\n    seis_data = []\n    source_x_locations = [0, 17, 34, 52, 69] # Using 0-based indices now\n    \n    for sx_idx in source_x_locations:\n        coord['sx'] = sx_idx * dx\n        \n        if method == 'abc24':\n            seis = a2d_mod_abc24(vel, nbc, dx, nt, dt, s, coord, isFS)\n        else:\n            raise ValueError(f\"Invalid method: {method}\")\n\n        seis_data.append(seis)\n        \n    return np.stack(seis_data, axis=0)\n    \n    \ndef plot_seis(seis):\n    \"\"\"\n    seis : (5, 1000, 70)\n    \"\"\"\n    fig,ax=plt.subplots(1,5,figsize=(20,5))\n    ax[0].imshow(seis[0, :, :],extent=[0,70,1000,0],aspect='auto',cmap='gray',vmin=-0.5,vmax=0.5)\n    ax[1].imshow(seis[1, :, :],extent=[0,70,1000,0],aspect='auto',cmap='gray',vmin=-0.5,vmax=0.5)\n    ax[2].imshow(seis[2, :, :],extent=[0,70,1000,0],aspect='auto',cmap='gray',vmin=-0.5,vmax=0.5)\n    ax[3].imshow(seis[3, :, :],extent=[0,70,1000,0],aspect='auto',cmap='gray',vmin=-0.5,vmax=0.5)\n    ax[4].imshow(seis[4, :, :],extent=[0,70,1000,0],aspect='auto',cmap='gray',vmin=-0.5,vmax=0.5)\n    for axis in ax:\n        axis.set_xticks(range(0, 70, 10))\n        axis.set_xticklabels(range(0, 700, 100))\n        axis.set_yticks(range(0, 2000, 1000))\n        axis.set_yticklabels(range(0, 2,1))\n        axis.set_ylabel('Time (s)', fontsize=12)\n        axis.set_xlabel('Offset (m)', fontsize=12)\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-07T07:58:39.067035Z","iopub.execute_input":"2025-06-07T07:58:39.067592Z","iopub.status.idle":"2025-06-07T07:58:39.096066Z","shell.execute_reply.started":"2025-06-07T07:58:39.067562Z","shell.execute_reply":"2025-06-07T07:58:39.094492Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"------\n## **Comparison**\nWhen I enter the velocity map provided by OpenFWI into the converted python code, I will check if Seismic data is generated well.  \nThe generated Seismic data and the Seismic data provided by OpenFWI should be almost identical.","metadata":{}},{"cell_type":"code","source":"seis_files = [\n    \"/kaggle/input/waveform-inversion/train_samples/FlatVel_A/data/data1.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/FlatVel_B/data/data1.npy\",\n    \n    \"/kaggle/input/waveform-inversion/train_samples/CurveVel_A/data/data1.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/CurveVel_B/data/data1.npy\",\n\n    \"/kaggle/input/waveform-inversion/train_samples/FlatFault_A/seis2_1_0.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/FlatFault_B/seis6_1_0.npy\",\n    \n    \"/kaggle/input/waveform-inversion/train_samples/CurveFault_A/seis2_1_0.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/CurveFault_B/seis6_1_0.npy\",\n\n    \"/kaggle/input/waveform-inversion/train_samples/Style_A/data/data1.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/Style_B/data/data1.npy\",\n    \n]\nvel_files = [\n    \"/kaggle/input/waveform-inversion/train_samples/FlatVel_A/model/model1.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/FlatVel_B/model/model1.npy\",\n    \n    \"/kaggle/input/waveform-inversion/train_samples/CurveVel_A/model/model1.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/CurveVel_B/model/model1.npy\",\n\n    \"/kaggle/input/waveform-inversion/train_samples/FlatFault_A/vel2_1_0.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/FlatFault_B/vel6_1_0.npy\",\n    \n    \"/kaggle/input/waveform-inversion/train_samples/CurveFault_A/vel2_1_0.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/CurveFault_B/vel6_1_0.npy\",\n\n    \"/kaggle/input/waveform-inversion/train_samples/Style_A/model/model1.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/Style_B/model/model1.npy\",\n    \n]\ntypes = [\n    'FlatVel_A', 'FlatVel_B', 'CurveVel_A', 'CurveVel_B', 'FlatFault_A', 'FlatFault_B', 'CurveFault_A', 'CurveFault_B', 'Style_A', 'Style_B'\n]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-07T07:58:39.098136Z","iopub.execute_input":"2025-06-07T07:58:39.098428Z","iopub.status.idle":"2025-06-07T07:58:39.126934Z","shell.execute_reply.started":"2025-06-07T07:58:39.098406Z","shell.execute_reply":"2025-06-07T07:58:39.125572Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nfor i, tp in enumerate(types):\n    seis_data = np.load(seis_files[i])\n    vel_data = np.load(vel_files[i])\n    seis_data_sim = vel_to_seis(vel_data[0][0])\n    print(f\"TYPE : {tp}\")\n    print(\"\\n< VELOCITY MAP >\")\n    fig, ax = plt.subplots(1, 1, figsize=(5, 2.5))\n    ax.imshow(vel_data[0, 0])\n    plt.show()\n\n    print(f\"\\n< FORWARD SIMULATION (from velocity map) : shape={seis_data_sim.shape} >\")\n    plot_seis(seis_data_sim[:, 1:, :])\n    \n    print(\"\\n< INPUT (origin) >\")\n    plot_seis(seis_data[0])\n    \n    print(\"\\n< ERROR >\")\n    errors = []\n    for j in range(5):\n        error = np.mean(np.abs(seis_data[0, j, :, :] - seis_data_sim[j, 1:, :])) # MAE\n        errors.append(error)\n        print(f\"Receiver{j+1} Error : {error:.6f}\")\n    print(f\"Mean Error : {np.mean(errors):.6f}\")\n    print(\"#########################################\")\n    print()\n\n    # if i==2:\n    #     break\n    # break","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-07T07:58:39.128457Z","iopub.execute_input":"2025-06-07T07:58:39.129009Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null}]}