Merge pull request #17 from AmirMardan/16-error-occurs-in-the-simple-…

…fwi-example fixing the device selection and a bug in `tauxz` buffer
AmirMardan · Nov 5, 2024 · 93e6b88 · 93e6b88
2 parents c04255f + 894cdb8
commit 93e6b88
Show file tree

Hide file tree

Showing 25 changed files with 2,158 additions and 258 deletions.
diff --git a/.github/workflows/pr.yml b/.github/workflows/pr.yml
@@ -0,0 +1,41 @@
+# This workflow will install Python dependencies, run tests and lint with a variety of Python versions
+# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python
+
+name: Python package
+
+on:
+  push:
+    branches: [ "dev" ]
+  pull_request:
+    branches: [ "main" ]
+
+jobs:
+  build:
+
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-latest, macos-latest, windows-latest]
+        python-version: ["3.8", "3.9", "3.10"]
+
+    steps:
+    - uses: actions/checkout@v3
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v3
+      with:
+        python-version: ${{ matrix.python-version }}
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        python -m pip install flake8 pytest
+        pip install -r requirements.txt
+    - name: Lint with flake8
+      run: |
+        # stop the build if there are Python syntax errors or undefined names
+        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
+        # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
+        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+    - name: Test with pytest
+      run: |
+        pytest
diff --git a/.gitignore b/.gitignore
@@ -0,0 +1,2 @@
+*.pyc
+.DS_Store
diff --git a/README.md b/README.md
@@ -7,6 +7,7 @@
 <img src="https://visitor-badge.laobi.icu/badge?page_id=pyfwi_unique12.pyfwi_unique12"/>
 
 This repository contains Python package for elastic seismic full-waveform inversion (FWI) and time-lapse FWI.
+Documentation of PyFWI is available [here](https://pyfwi.readthedocs.io/en/latest/index.html).
 
 
 ## Installation

diff --git a/docs/CNAME b/docs/CNAME
@@ -0,0 +1 @@
+pyfwi.readthedocs.io
diff --git a/docs/conf.py b/docs/conf.py
@@ -20,7 +20,7 @@
 project = 'PyFWI'
 copyright = '2022, Amir Mardan'
 author = 'Amir Mardan'
-release = '0.1.9'
+release = '0.1.10'
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration

diff --git a/docs/index.rst b/docs/index.rst
@@ -9,6 +9,13 @@ Welcome to PyFWI's documentation!
 **PyFWI** is an open source Python package to perform seismic modeling and full-waveform inversion (FWI) in elastic media. 
 This package is implemented in time domain and coded using GPU programming (PyOpenCL) to accelerate the computation.
 
+If you have any questions about PyFWI, please use the following this `link <https://github.com/AmirMardan/PyFWI/discussions>`_.
+
+If you have any technical questions about FWI, TL-FWI, or PyFWI, please visit `my personal website <https://amirmardan.github.io/tlfwi.html/>`_.
+All my publications are available there. I will be happy to assist if you contact me via email.
+
+For bugs, developments, and errors, please use issues in the GitHub repository available `here <https://github.com/AmirMardan/PyFWI/issues>`_ to ask your questions.
+
 .. toctree::
    :maxdepth: 1
    :caption: Getting started
@@ -62,4 +69,4 @@ Citing PyFWI
   year = {2023},
   publisher = {Elsevier},
   doi = {10.1016/j.softx.2023.101384}
-}
+   }
diff --git a/docs/sub_doc/example.rst b/docs/sub_doc/example.rst
@@ -73,8 +73,7 @@ Then we need to create an input dictionary as follow
         '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
+        'pml_dir': 2,  # type of boundary layer 
     }
     
     seisout = 0 # Type of output 0: Pressure

diff --git a/docs/sub_doc/fwi_example.rst b/docs/sub_doc/fwi_example.rst
@@ -56,7 +56,6 @@ Then we need to create an input dictionary as follow
         '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
@@ -193,149 +192,18 @@ these parameters as
    k_0 = 1
    k_end = 3
 
+Let's call the FWI object,
+
 .. code:: ipython3
 
     m_est, _ = fwi(m0, method="lbfgs", 
                      freqs=[25, 45], iter=[2, 2], 
                      n_params=1, k_0=1, k_end=2)
 
-
-.. parsed-literal::
-
-    Parameter number 1 to 1
-    2500.0 2500.0
-     for f= 25: rms is: 0.0003187612455803901 with rms_reg: 0, and rms_data: 0.0003187612455803901, rms_mp: 0.0, rms_model_relation: 0
-    Parameter number 1 to 1
-    RUNNING THE L-BFGS-B CODE
-    
-               * * *
-    
-    Machine precision = 2.220D-16
-     N =        10000     M =           10
-    
-    At X0         0 variables are exactly at the bounds
-    
-    At iterate    0    f=  3.18761D-04    |proj g|=  7.19684D-10
-    
-               * * *
-    
-    Tit   = total number of iterations
-    Tnf   = total number of function evaluations
-    Tnint = total number of segments explored during Cauchy searches
-    Skip  = number of BFGS updates skipped
-    Nact  = number of active bounds at final generalized Cauchy point
-    Projg = norm of the final projected gradient
-    F     = final function value
-    
-               * * *
-    
-       N    Tit     Tnf  Tnint  Skip  Nact     Projg        F
-    10000      0      1      0     0     0   7.197D-10   3.188D-04
-      F =   3.1876124558039010E-004
-    
-    CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL            
-
-
-.. parsed-literal::
-
-     This problem is unconstrained.
-
-
-.. parsed-literal::
-
-    2500.0 2500.0
-     for f= 45: rms is: 0.004415073432028294 with rms_reg: 0, and rms_data: 0.004415073432028294, rms_mp: 0.0, rms_model_relation: 0
-    RUNNING THE L-BFGS-B CODE
-    
-               * * *
-    
-    Machine precision = 2.220D-16
-     N =        10000     M =           10
-    
-    At X0         0 variables are exactly at the bounds
-    
-    At iterate    0    f=  4.41507D-03    |proj g|=  2.39370D-08
-    
-    
-    ITERATION     1
-    
-    ---------------- CAUCHY entered-------------------
-     There are            0   breakpoints 
-    
-     GCP found in this segment
-    Piece      1 --f1, f2 at start point  -6.4549D-13  6.4549D-13
-    Distance to the stationary point =   1.0000D+00
-    
-    ---------------- exit CAUCHY----------------------
-    
-           10000  variables are free at GCP            1
-
-
-.. parsed-literal::
-
-     This problem is unconstrained.
-
-
-.. parsed-literal::
-
-    2499.9827111341883 2500.0297937900996
-     for f= 45: rms is: 0.004414360038936138 with rms_reg: 0, and rms_data: 0.004414360038936138, rms_mp: 0.0, rms_model_relation: 0
-    2499.9135556709416 2500.1489689504983
-     for f= 45: rms is: 0.004411510191857815 with rms_reg: 0, and rms_data: 0.004411510191857815, rms_mp: 0.0, rms_model_relation: 0
-    2499.636933817955 2500.6256695920933
-     for f= 45: rms is: 0.004400133620947599 with rms_reg: 0, and rms_data: 0.004400133620947599, rms_mp: 0.0, rms_model_relation: 0
-    2498.530446406008 2502.5324721584725
-     for f= 45: rms is: 0.00435509392991662 with rms_reg: 0, and rms_data: 0.00435509392991662, rms_mp: 0.0, rms_model_relation: 0
-    2494.1044967582216 2510.1596824239887
-     for f= 45: rms is: 0.004182371310889721 with rms_reg: 0, and rms_data: 0.004182371310889721, rms_mp: 0.0, rms_model_relation: 0
-    2476.400698167074 2540.668523486055
-     for f= 45: rms is: 0.003607903141528368 with rms_reg: 0, and rms_data: 0.003607903141528368, rms_mp: 0.0, rms_model_relation: 0
-     LINE SEARCH           5  times; norm of step =    1365.0000000000000     
-    
-    At iterate    1    f=  3.60790D-03    |proj g|=  1.91296D-08
-    
-    
-    ITERATION     2
-    
-    ----------------SUBSM entered-----------------
-    
-    
-    ----------------exit SUBSM --------------------
-    
-    2356.1934399025804 2683.887118020453
-     for f= 45: rms is: 0.0026031285524368286 with rms_reg: 0, and rms_data: 0.0026031285524368286, rms_mp: 0.0, rms_model_relation: 0
-     LINE SEARCH           0  times; norm of step =    3577.6768283163833     
-    
-    At iterate    2    f=  2.60313D-03    |proj g|=  1.26216D-08
-    
-               * * *
-    
-    Tit   = total number of iterations
-    Tnf   = total number of function evaluations
-    Tnint = total number of segments explored during Cauchy searches
-    Skip  = number of BFGS updates skipped
-    Nact  = number of active bounds at final generalized Cauchy point
-    Projg = norm of the final projected gradient
-    F     = final function value
-    
-               * * *
-    
-       N    Tit     Tnf  Tnint  Skip  Nact     Projg        F
-    10000      2      8      1     0     0   1.262D-08   2.603D-03
-      F =   2.6031285524368286E-003
-    
-    STOP: TOTAL NO. of ITERATIONS REACHED LIMIT                 
-
-
-Here is the estimated model
-
 .. code:: ipython3
 
     # Time to plot the results
     fig = splt.earth_model(m_est, ['vp'], cmap='jet')
 
-
-
-
 .. image:: fwi_example_files/fwi_example_21_0.png
 
diff --git a/docs/sub_doc/fwi_example_files/fwi_example_15_1.png b/docs/sub_doc/fwi_example_files/fwi_example_15_1.png
diff --git a/docs/sub_doc/grad_pytorch.rst b/docs/sub_doc/grad_pytorch.rst
@@ -66,8 +66,7 @@ Then we need to create an input dictionary as follow
         '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
-        'seimogram_shape': '3d',
+        'seimogram_shape': '3d' 
     }
     
     seisout = 0 # Type of output 0: Pressure

diff --git a/example/fwi_example.ipynb → example/fwi_example_crosswell.ipynb b/example/fwi_example.ipynb → example/fwi_example_crosswell.ipynb
diff --git a/example/fwi_example_surface.ipynb b/example/fwi_example_surface.ipynb
diff --git a/example/fwi_example_surface.py b/example/fwi_example_surface.py
@@ -0,0 +1,95 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+
+import sys
+sys.path.append(os.path.abspath('./src/'))
+
+import PyFWI.wave_propagation as wave
+import PyFWI.acquisition as acq
+import PyFWI.seiplot as splt
+import PyFWI.model_dataset as md
+from PyFWI.fwi import FWI
+
+Model = md.ModelGenerator('louboutin')
+model = Model()
+
+im = splt.earth_model(model, cmap='coolwarm')
+
+model_shape = model[[*model][0]].shape
+#%% ================== Parameters and Geometry ==================
+
+inpa = {
+    'ns': 4,  # 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.6,  # 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
+}
+
+seisout = 0 # Type of output 0: Pressure
+
+inpa['rec_dis'] =  1 * inpa['dh']  # Define the receivers' distance
+
+
+offsetx = inpa['dh'] * model_shape[1]
+depth = inpa['dh'] * model_shape[0]
+
+src_loc, rec_loc, n_surface_rec, n_well_rec = acq.acq_parameters(inpa['ns'], 
+                                                                 inpa['rec_dis'], 
+                                                                 offsetx,
+                                                                 depth,
+                                                                 inpa['dh'], 
+                                                                 inpa['sdo'], 
+                                                                 acq_type=inpa['acq_type'])        
+rec_loc[:, 1] -= 2 * inpa['dh']
+
+# Create the source
+src = acq.Source(src_loc, inpa['dh'], inpa['dt'])
+src.Ricker(inpa['fdom'])
+
+#%% ================== Forward Modelling ==================
+# Create the wave object
+W = wave.WavePropagator(inpa, src, rec_loc, model_shape,
+                        n_well_rec=n_well_rec,
+                        components=seisout, chpr=0)
+
+# Call the forward modelling 
+d_obs = W.forward_modeling(model, show=False)  # show=True can show the propagation of the wave
+
+plt.imshow(d_obs["taux"], cmap='gray', 
+           aspect="auto", extent=[0, d_obs["taux"].shape[1], inpa['t'], 0])
+
+#%% ================== Initial model ==================
+m0 = Model(smoothing=1)
+m0['vs'] *= 0.0
+m0['rho'] = np.ones_like(model['rho'])
+
+fig = splt.earth_model(m0, ['vp'], cmap='coolwarm')
+
+fig.axes[0].plot(src_loc[:,0]//inpa["dh"], 
+                 src_loc[:,1]//inpa["dh"], "rv", markersize=5)
+
+fig.axes[0].plot(rec_loc[:,0]//inpa["dh"], 
+                 rec_loc[:,1]//inpa["dh"], "b*", markersize=3)
+
+
+#%% ================== FWI ==================
+inpa['energy_balancing'] = True
+
+fwi = FWI(d_obs, inpa, src, rec_loc, model_shape, 
+          components=seisout, chpr=20, n_well_rec=n_well_rec)
+
+m_est, _ = fwi(m0, method="lbfgs", 
+                 freqs=[25, 45], iter=[2, 2], 
+                 n_params=1, k_0=1, k_end=2)
+
+fig = splt.earth_model(m_est, ['vp'], cmap='jet')
+
+plt.show()
+a = 1
diff --git a/example/grad_pytorch.ipynb b/example/grad_pytorch.ipynb
diff --git a/example/gradient_example.ipynb b/example/gradient_example.ipynb