add documentation for fwi

docs/index.rst

@@ -14,6 +14,7 @@ This package is implemented in time domain and coded using GPU programming (PyOp
    :caption: Getting started
 
    sub_doc/example
+   sub_doc/fwi_example
 
 .. toctree::
    :maxdepth: 1
@@ -52,14 +53,13 @@ Citing PyFWI
 
 ::
 
-   @software{PyFWI,
-   author       = {Mardan Amir and
-                  Bernard Giroux and
-                  Gabriel Fabien-Ouellet},
-   title        = {{PyFWI}: A {Python} package for full-waveform inversion and reservoir monitoring},
-   month        = Jan,
-   year         = 2022,
-   publisher    = {Zenodo},
-   doi          = {10.5281/zenodo.5813637},
-   url          = {https://doi.org/10.5281/zenodo.5813637}
-   }
+   @article{mardan2023pyfwi,
+  title = {PyFWI: {A Python} package for full-waveform inversion and reservoir monitoring},
+  author = {Mardan, Amir and Giroux, Bernard and Fabien-Ouellet, Gabriel},
+  journal = {SoftwareX},
+  volume = {22},
+  pages = {101384},
+  year = {2023},
+  publisher = {Elsevier},
+  doi = {10.1016/j.softx.2023.101384}
+}

docs/sub_doc/example.rst

@@ -0,0 +1,242 @@
+Installation
+============
+
+PyFWI can be installed using ``pip`` as
+
+.. code:: console
+
+
+       (.venv) $ pip install PyFWI
+
+on macOS or
+
+.. code:: console
+
+
+       (.venv) $ py -m pip install PyFWI
+
+on Windows.
+
+Simple Gradient Computation
+============================
+
+In this section we see some applications of PyFWI. First, forward
+modeling is shown and then we estimate gradient of cost funtion with
+respect to :math:`V_P`.
+
+**1. Forward modeling**
+
+In this simple example, we use PyFWI to do forward modeling. So, we need
+to first import the following packages amd modulus.
+
+.. code:: ipython3
+
+    import matplotlib.pyplot as plt
+    import numpy as np
+    
+    import PyFWI.wave_propagation as wave
+    import PyFWI.acquisition as acq
+    import PyFWI.seiplot as splt
+    import PyFWI.model_dataset as md
+    import PyFWI.fwi_tools as tools
+    import PyFWI.processing as process
+    
+
+
+A simple model can be created by using ``model_dataset`` module as
+
+.. code:: ipython3
+
+    Model = md.ModelGenerator('louboutin')
+    model = Model()
+    
+    im = splt.earth_model(model, cmap='coolwarm')
+
+
+
+.. image:: example_files/example_7_0.png
+
+
+Then we need to create an input dictionary as follow
+
+.. code:: ipython3
+
+    model_shape = model[[*model][0]].shape
+    
+    inpa = {
+        'ns': 1,  # Number of sources
+        'sdo': 4,  # Order of FD
+        'fdom': 15,  # Central frequency of source
+        'dh': 7,  # Spatial sampling rate
+        'dt': 0.004,  # Temporal sampling rate
+        'acq_type': 1,  # Type of acquisition (0: crosswell, 1: surface, 2: both)
+        't': 0.8,  # Length of operation
+        'npml': 20,  # Number of PML 
+        'pmlR': 1e-5,  # Coefficient for PML (No need to change)
+        'pml_dir': 2,  # type of boundary layer
+        'device': 1, # The device to run the program. Usually 0: CPU 1: GPU
+    }
+    
+    seisout = 0 # Type of output 0: Pressure
+    
+    inpa['rec_dis'] =  1 * inpa['dh']  # Define the receivers' distance
+
+
+Now, we obtain the location of sources and receivers based on specified
+parameters.
+
+.. code:: ipython3
+
+    offsetx = inpa['dh'] * model_shape[1]
+    depth = inpa['dh'] * model_shape[0]
+    
+    src_loc, rec_loc = acq.surface_seismic(inpa['ns'], inpa['rec_dis'], offsetx,
+                                                          inpa['dh'], inpa['sdo'])        
+    src_loc[:, 1] -= 5 * inpa['dh']
+    
+    # Create the source
+    src = acq.Source(src_loc, inpa['dh'], inpa['dt'])
+    src.Ricker(inpa['fdom'])
+
+
+Finally, we can have the forward modelling as
+
+.. code:: ipython3
+
+    # Create the wave object
+    W = wave.WavePropagator(inpa, src, rec_loc, model_shape, components=seisout, chpr=20)
+    
+    # Call the forward modelling 
+    d_obs = W.forward_modeling(model, show=False)  # show=True can show the propagation of the wave
+
+To compute the gradient using the adjoint-state method, we need to save
+the wavefield during the forward wave propagation. This must be done for
+the wavefield obtained from estimated model. For example, the wavefield
+at four time steps are presented here in addition to a shot gather.
+
+.. code:: ipython3
+
+    fig = plt.figure(figsize=(8, 4))
+    
+    count = 1
+    
+    ax = fig.add_subplot(122)
+    ax = splt.seismic_section(ax, d_obs['taux'], t_axis=np.linspace(0, inpa['t'], int(1 + inpa['t'] // inpa['dt'])))
+    
+    ax_loc = [1, 2, 5, 6]
+    snapshots = [40, 80, 130, 180]
+    
+    for i in range(len(snapshots)):
+        ax = fig.add_subplot(2, 4, ax_loc[i])
+        ax.imshow(W.W['taux'][:, :, 0, snapshots[i]], cmap='coolwarm')
+        
+        ax.axis('off')
+        count += 1
+    fig.suptitle("Wave propagation and a shot gather", fontweight='bold');
+
+
+
+
+.. image:: example_files/example_15_0.png
+
+
+**2. Gradient**
+
+To compute the gradient, we need the observed data and an initial model.
+So, first we obtain the observed data using more sources.
+
+**Note:** For better visualization and avoiding crosstalk, I compute the
+gradient in acoustic media.
+
+.. code:: ipython3
+
+    # Making medium acoustic
+    model['vs'] *= 0.0
+    model['rho'] = np.ones_like(model['rho'])
+    
+    # Increasing number of sources
+    inpa['ns'] = 5
+    
+    src_loc, rec_loc = acq.surface_seismic(inpa['ns'], inpa['rec_dis'], offsetx,
+                                                          inpa['dh'], inpa['sdo'])        
+    src_loc[:, 1] -= 5 * inpa['dh']
+    
+    # Create the source
+    src = acq.Source(src_loc, inpa['dh'], inpa['dt'])
+    src.Ricker(inpa['fdom'])
+    
+    # Create the wave object
+    W = wave.WavePropagator(inpa, src, rec_loc, model_shape, components=seisout, chpr=20)
+    
+    # Call the forward modelling 
+    db_obs = W.forward_modeling(model, show=False)  # show=True can show the propagation of the wave
+    
+    # preparing data amd applying gain if required
+    db_obs = process.prepare_residual(db_obs, 1)
+
+Then we create the initial model.
+
+.. code:: ipython3
+
+    m0 = Model(smoothing=1)
+    m0['vs'] *= 0.0
+    m0['rho'] = np.ones_like(model['rho'])
+    
+    im = splt.earth_model(m0, ['vp'], cmap='coolwarm')
+
+
+
+.. image:: example_files/example_20_0.png
+
+
+And we simulate the wave propagation to obtain estimated data. For
+computing the gradient, we can smooth the gradient and scale it by
+defining ``g_smooth`` and ``energy_balancing``.
+
+.. code:: ipython3
+
+    inpa['energy_balancing'] = True
+
+We save the wavefield at 20% of the time steps (``chpr = 20``) to be
+used for gradient calculation. The value of wavefield is accessible
+using the attribute ``W`` which is a dictionary for :math:`V_x`,
+:math:`V_z`, :math:`\tau_x`, :math:`\tau_z`, and :math:`\tau_{xz}` as
+``vx``, ``vz``, ``taux``, ``tauz``, and ``tauxz``. Each parameter is a
+4D tensor. For example, we can have access to the last time step of
+:math:`\tau_x` for the first shot as ``W.W['taux'][:, :, 0, -1]``.
+
+.. code:: ipython3
+
+    Lam = wave.WavePropagator(inpa, src, rec_loc, model_shape,
+                              chpr=20, components=seisout)
+    
+    d_est = Lam.forward_modeling(m0, False)
+    d_est = process.prepare_residual(d_est, 1)
+
+
+Now, we define the cost function and obtaine the residuals for
+adjoint-state method.
+
+.. code:: ipython3
+
+    CF = tools.CostFunction('l2')
+    rms, adj_src = tools.cost_seismic(d_est, db_obs, fun=CF)
+    # print(rms)
+
+Using the adjoint source, we can estimate the gradient as
+
+.. code:: ipython3
+
+    grad = Lam.gradient(adj_src, show=False)
+
+.. code:: ipython3
+
+    # Time to plot the results
+    splt.earth_model(grad, ['vp'], cmap='jet');
+
+
+
+
+.. image:: example_files/example_29_0.png
+
+