Updated code to adopt Rommel's notation

groupeLIAMG · Feb 7, 2024 · 7deeada · 7deeada
1 parent 4f6a017
commit 7deeada
Show file tree

Hide file tree

Showing 34 changed files with 2,263 additions and 446 deletions.
diff --git a/examples/example2.ipynb b/examples/example2.ipynb
diff --git a/examples/example2.vtu b/examples/example2.vtu
diff --git a/examples/example4.ipynb b/examples/example4.ipynb
diff --git a/examples/example5.ipynb b/examples/example5.ipynb
diff --git a/examples/figs/example2_model.pdf b/examples/figs/example2_model.pdf
diff --git a/examples/figs/example2_rays.pdf b/examples/figs/example2_rays.pdf
diff --git a/examples/figs/example2_rays2.pdf b/examples/figs/example2_rays2.pdf
diff --git a/examples/figs/example2_tt.pdf b/examples/figs/example2_tt.pdf
diff --git a/src/ttcrpy/rgrid.pxd b/src/ttcrpy/rgrid.pxd
@@ -158,8 +158,8 @@ cdef extern from "Grid2D.h" namespace "ttcr" nogil:
         void setDelta(vector[T1]&) except +
         void setEpsilon(vector[T1]&) except +
         void setGamma(vector[T1]&) except +
-        void setR2(vector[T1]&) except +
-        void setR4(vector[T1]&) except +
+        void setS2(vector[T1]&) except +
+        void setS4(vector[T1]&) except +
         T1 computeSlowness(S&) except +
         void getTT(vector[T1]& tt, size_t threadNo) except +
         void raytrace(vector[S]& Tx,

diff --git a/src/ttcrpy/rgrid.pyx b/src/ttcrpy/rgrid.pyx
@@ -1832,7 +1832,7 @@ cdef class Grid2d:
             - 'tilted_elliptical' : tilted elliptical anisotropy
             - 'vti_psv' : vertical transverse isotropy, P and SV waves
             - 'vti_sh' : vertical transverse isotropy, SH waves
-            - 'weakly_anelliptical' : Weakly-Anelliptical
+            - 'weakly_anelliptical' : Weakly-Anelliptical formulation of B. Rommel
     eps : double
         convergence criterion (FSM) (default is 1e-15)
     maxit : int
@@ -2499,11 +2499,11 @@ cdef class Grid2d:
             raise ValueError('v must be 1D or 2D ndarray')
         self.grid.setGamma(data)
 
-    def set_r2(self, v):
+    def set_s2(self, v):
         """
-        set_r2(g)
+        set_s2(g)
 
-        Assign weakly anelliptical parameter r2
+        Assign weakly anelliptical parameter s2
 
         Parameters
         ----------
@@ -2533,13 +2533,13 @@ cdef class Grid2d:
                 data.push_back(tmp[i])
         else:
             raise ValueError('v must be 1D or 2D ndarray')
-        self.grid.setR2(data)
+        self.grid.setS2(data)
 
-    def set_r4(self, v):
+    def set_s4(self, v):
         """
-        set_r4(g)
+        set_s4(g)
 
-        Assign weakly anelliptical parameter r4
+        Assign weakly anelliptical parameter s4
 
         Parameters
         ----------
@@ -2569,7 +2569,7 @@ cdef class Grid2d:
                 data.push_back(tmp[i])
         else:
             raise ValueError('v must be 1D or 2D ndarray')
-        self.grid.setR4(data)
+        self.grid.setS4(data)
 
     def compute_D(self, np.ndarray[np.double_t, ndim=2] coord):
         """

diff --git a/src/ttcrpy/tmesh.pxd b/src/ttcrpy/tmesh.pxd
@@ -165,8 +165,8 @@ cdef extern from "Grid2D.h" namespace "ttcr" nogil:
         void getSlowness(vector[T1]&) except +
         void setXi(vector[T1]&) except +
         void setTiltAngle(vector[T1]&) except +
-        void setR2(vector[T1]&) except +
-        void setR4(vector[T1]&) except +
+        void setS2(vector[T1]&) except +
+        void setS4(vector[T1]&) except +
         void getTT(vector[T1]& tt, size_t threadNo) except +
         void raytrace(vector[S]& Tx,
                       vector[T1]& t0,

diff --git a/src/ttcrpy/tmesh.pyx b/src/ttcrpy/tmesh.pyx
@@ -1216,7 +1216,7 @@ cdef class Mesh2d:
             - 'iso' : isotropic medium
             - 'elliptical' : elliptical anisotropy
             - 'tilted_elliptical' : tilted elliptical anisotropy
-            - 'weakly_anelliptical' : Weakly-Anelliptical
+            - 'weakly_anelliptical' : Weakly-Anelliptical formulation of B. Rommel
     eps : double
         convergence criterion (FSM) (default is 1e-15)
     maxit : int
@@ -1597,49 +1597,49 @@ cdef class Mesh2d:
             var.push_back(theta[i])
         self.grid.setTiltAngle(var)
 
-    def set_r2(self, r2):
+    def set_s2(self, s2):
         """
-        set_r2(r2)
+        set_s2(s2)
 
-        Assign phase-velocity parameter r2 to grid
+        Assign energy-velocity parameter s2 to grid
 
         Parameters
         ----------
-        r2 : np ndarray, shape (nparams, )
+        s2 : np ndarray, shape (nparams, )
         """
-        if r2.size != self.nparams:
-            raise ValueError('r2 vector has wrong size')
+        if s2.size != self.nparams:
+            raise ValueError('s2 vector has wrong size')
 
-        if not r2.flags['C_CONTIGUOUS']:
-                r2 = np.ascontiguousarray(r2)
+        if not s2.flags['C_CONTIGUOUS']:
+                s2 = np.ascontiguousarray(s2)
 
         cdef vector[double] var
         cdef int i
-        for i in range(r2.size):
-            var.push_back(r2[i])
-        self.grid.setR2(var)
+        for i in range(s2.size):
+            var.push_back(s2[i])
+        self.grid.setS2(var)
 
-    def set_r4(self, r4):
+    def set_s4(self, s4):
         """
-        set_r4(r4)
+        set_s4(s4)
 
-        Assign phase-velocity parameter r4 to grid
+        Assign energy-velocity parameter s4 to grid
 
         Parameters
         ----------
-        r4 : np ndarray, shape (nparams, )
+        s4 : np ndarray, shape (nparams, )
         """
-        if r4.size != self.nparams:
-            raise ValueError('r4 vector has wrong size')
+        if s4.size != self.nparams:
+            raise ValueError('s4 vector has wrong size')
 
-        if not r4.flags['C_CONTIGUOUS']:
-                r4 = np.ascontiguousarray(r4)
+        if not s4.flags['C_CONTIGUOUS']:
+                s4 = np.ascontiguousarray(s4)
 
         cdef vector[double] var
         cdef int i
-        for i in range(r4.size):
-            var.push_back(r4[i])
-        self.grid.setR4(var)
+        for i in range(s4.size):
+            var.push_back(s4[i])
+        self.grid.setS4(var)
 
     def raytrace(self, source, rcv, slowness=None, thread_no=None,
                  aggregate_src=False, compute_L=False, return_rays=False):

diff --git a/tests/files/Grid2Drcsp_tt_grid_elliptical.vtr b/tests/files/Grid2Drcsp_tt_grid_elliptical.vtr
@@ -10,13 +10,13 @@
     <CellData>
     </CellData>
     <Coordinates>
-      <DataArray type="Float64" Name="Array 0x10620e020" format="binary" RangeMin="0" RangeMax="6000">
+      <DataArray type="Float64" Name="Array 0x1132195d0" format="binary" RangeMin="0" RangeMax="6000">
         AQAAAACAAAAoAwAAFwEAAA==eJxNzi1MQmEUxvEbCQYiwUAkGAgG4g0GEiMYSI7AHIndOeacwd0x5phzehVUBJUroCL4cf2EQCAaDEQD4UaDgUgwGJ5/8C2/ved9zzmPYfw/CVOuSHuVe06O1qivS3OT9y3p2/wrSLfI/22ZLtG3I8O79O/Rv4+OHB8w71B6ZeZWpHPE/GNpnbCnKpOn7KvJaJ29ZzJ4zv4LOUW/IbMud0xdkgvjTfJhrEVOjLTJi6ErcmPgmvw4Q+tGfmO6I78weSs/0OzKAUZ7soPhO1nF4L0sofEgN3CK2Ufp44JHHdvo4/yTTGEFxzj3LONYxBH+YuxF5tHDH4y8ygy6OMHQm1xGBz8x8C6X0MYhznCxL62++QcopYc0
       </DataArray>
-      <DataArray type="Float64" Name="Array 0x106207550" format="binary" RangeMin="0" RangeMax="0">
+      <DataArray type="Float64" Name="Array 0x11321a580" format="binary" RangeMin="0" RangeMax="0">
         AQAAAACAAAAIAAAACwAAAA==eJxjYIAAAAAIAAE=
       </DataArray>
-      <DataArray type="Float64" Name="Array 0x106208500" format="binary" RangeMin="0" RangeMax="3000">
+      <DataArray type="Float64" Name="Array 0x11321a6c0" format="binary" RangeMin="0" RangeMax="3000">
         AQAAAACAAACYAQAAqQAAAA==eJxNzikOwkAYhuFKJLISiUQiR6JIJYqgCIo0CASCNIQQQggtW9lpS1ksR+AIHKFHQCIR3ysY82Rm/uWzrP9TNbIuvSb3tnx1eO9K0+O/LzOPuoGMh9SPZGNM30QWpvTP6PcxkO858xbyuWTuSgZr5ofS3bBnK50d+/aydGDvUeZP7D/LD2aRbMXcsZaQCysX8mE5JScWr+RF+0ZuzN3Jj190H+YHwvo/dQ==
       </DataArray>
     </Coordinates>

diff --git a/tests/files/Grid2Drcsp_tt_grid_weakly.vtr b/tests/files/Grid2Drcsp_tt_grid_weakly.vtr
diff --git a/tests/files/Grid2Ducsp_tt_grid_weakly.vtu b/tests/files/Grid2Ducsp_tt_grid_weakly.vtu
diff --git a/tests/files/elliptical_fine2d.vtr b/tests/files/elliptical_fine2d.vtr
@@ -13,13 +13,13 @@
       </DataArray>
     </CellData>
     <Coordinates>
-      <DataArray type="Float64" Name="Array 0x1306073f0" format="binary" RangeMin="0" RangeMax="6000">
+      <DataArray type="Float64" Name="Array 0x12d0b2f30" format="binary" RangeMin="0" RangeMax="6000">
         AQAAAACAAAAoAwAAFwEAAA==eJxNzi1MQmEUxvEbCQYiwUAkGAgG4g0GEiMYSI7AHIndOeacwd0x5phzehVUBJUroCL4cf2EQCAaDEQD4UaDgUgwGJ5/8C2/ved9zzmPYfw/CVOuSHuVe06O1qivS3OT9y3p2/wrSLfI/22ZLtG3I8O79O/Rv4+OHB8w71B6ZeZWpHPE/GNpnbCnKpOn7KvJaJ29ZzJ4zv4LOUW/IbMud0xdkgvjTfJhrEVOjLTJi6ErcmPgmvw4Q+tGfmO6I78weSs/0OzKAUZ7soPhO1nF4L0sofEgN3CK2Ufp44JHHdvo4/yTTGEFxzj3LONYxBH+YuxF5tHDH4y8ygy6OMHQm1xGBz8x8C6X0MYhznCxL62++QcopYc0
       </DataArray>
-      <DataArray type="Float64" Name="Array 0x130607530" format="binary" RangeMin="0" RangeMax="0">
+      <DataArray type="Float64" Name="Array 0x12d0b1780" format="binary" RangeMin="0" RangeMax="0">
         AQAAAACAAAAIAAAACwAAAA==eJxjYIAAAAAIAAE=
       </DataArray>
-      <DataArray type="Float64" Name="Array 0x130608990" format="binary" RangeMin="0" RangeMax="3000">
+      <DataArray type="Float64" Name="Array 0x12d09e3c0" format="binary" RangeMin="0" RangeMax="3000">
         AQAAAACAAACYAQAAqQAAAA==eJxNzikOwkAYhuFKJLISiUQiR6JIJYqgCIo0CASCNIQQQggtW9lpS1ksR+AIHKFHQCIR3ysY82Rm/uWzrP9TNbIuvSb3tnx1eO9K0+O/LzOPuoGMh9SPZGNM30QWpvTP6PcxkO858xbyuWTuSgZr5ofS3bBnK50d+/aydGDvUeZP7D/LD2aRbMXcsZaQCysX8mE5JScWr+RF+0ZuzN3Jj190H+YHwvo/dQ==
       </DataArray>
     </Coordinates>

diff --git a/tests/files/prepare_tests_aniso.py b/tests/files/prepare_tests_aniso.py
@@ -5,7 +5,7 @@
 import matplotlib.pyplot as plt
 import vtk
 
-import traveltime as wa_tt
+import ttwean
 
 # %%
 
@@ -241,15 +241,20 @@
 # writer.SetDataModeToBinary()
 # writer.Write()
 
-# %%
+# %% Weakly 
 
-VP0 = 1875.
-r2 = 0.1
-r4 = 0.15602005
+# Define "Dog Creek Shale" (Thomsen, 1986)
+vp0 = 1875.       # reference P-velocity
+vs0 = 826.        # reference S-Velocity
+delta = 0.100     # Thomsen's delta
+epsilon = 0.225   # Thomsen's epsilon
+
+vp00, vs00, sss = ttwean.TTWean.ttinit(vp0=vp0, vs0=vs0, delta=delta, epsilon=epsilon, wavetype='P')
 
-rrr = np.array([1., np.NAN, r2, np.NAN, r4])   # generic anisotropy parameters (here, for Dog Creek Shale)
+VP0 = 1875.
+s2 = sss[2]
+s4 = sss[4]
 
-scs = wa_tt.precompute(v00=VP0, rrr=rrr)                                 # precompute energy velocity and construct an interpolation polynomial
 
 x = np.arange(0.0, 6000.01, 60.0)
 z = np.arange(0.0, 3000.01, 60.0)
@@ -262,9 +267,14 @@
 xx = xx.flatten()
 zz = zz.flatten()
 
-t = wa_tt.traveltime(xxx=xx, zzz=zz, scs=scs)
+theta = np.arctan2(xx, zz)
+
+v = ttwean.TTWean.ttvel(vp0=vp0, vs0=vs0, sss=sss, wavetype='P', angle=theta)
+
+d = np.sqrt(xx*xx + zz*zz)
+t = d/v
 
-t[0] = 0.0
+#t[0] = 0.0
 
 
 # %%
@@ -287,19 +297,19 @@
 rgrid.SetZCoordinates(zCoords)
 
 v = vtk.vtkDoubleArray()
-r2_v = vtk.vtkDoubleArray()
-r4_v = vtk.vtkDoubleArray()
+s2_v = vtk.vtkDoubleArray()
+s4_v = vtk.vtkDoubleArray()
 
 for i in range(ncx*ncz):
     v.InsertNextValue(VP0)
-    r2_v.InsertNextValue(r2)
-    r4_v.InsertNextValue(r4)
+    s2_v.InsertNextValue(s2)
+    s4_v.InsertNextValue(s4)
 v.SetName("Velocity")
-r2_v.SetName("r2")
-r4_v.SetName("r4")
+s2_v.SetName("s2")
+s4_v.SetName("s4")
 rgrid.GetCellData().AddArray(v)
-rgrid.GetCellData().AddArray(r2_v)
-rgrid.GetCellData().AddArray(r4_v)
+rgrid.GetCellData().AddArray(s2_v)
+rgrid.GetCellData().AddArray(s4_v)
 
 writer = vtk.vtkXMLRectilinearGridWriter()
 writer.SetInputData(rgrid)
@@ -348,24 +358,24 @@
 
 triangle = vtk.vtkTriangle()
 v = vtk.vtkDoubleArray()
-r2_v = vtk.vtkDoubleArray()
-r4_v = vtk.vtkDoubleArray()
+s2_v = vtk.vtkDoubleArray()
+s4_v = vtk.vtkDoubleArray()
 v.SetName("Velocity")
-r2_v.SetName("r2")
-r4_v.SetName("r4")
+s2_v.SetName("s2")
+s4_v.SetName("s4")
 
 for n in range(len(tri)):
     triangle.GetPointIds().SetId(0, tri[n][0])
     triangle.GetPointIds().SetId(1, tri[n][1])
     triangle.GetPointIds().SetId(2, tri[n][2])
     ugrid.InsertNextCell( triangle.GetCellType(), triangle.GetPointIds() )
     v.InsertNextValue(VP0)
-    r2_v.InsertNextValue(r2)
-    r4_v.InsertNextValue(r4)
+    s2_v.InsertNextValue(s2)
+    s4_v.InsertNextValue(s4)
 
 ugrid.GetCellData().AddArray(v)
-ugrid.GetCellData().AddArray(r2_v)
-ugrid.GetCellData().AddArray(r4_v)
+ugrid.GetCellData().AddArray(s2_v)
+ugrid.GetCellData().AddArray(s4_v)
 
 writer = vtk.vtkXMLUnstructuredGridWriter()
 writer.SetInputData(ugrid)

diff --git a/tests/files/sol_analytique_elliptical_2d_tt.vtr b/tests/files/sol_analytique_elliptical_2d_tt.vtr
@@ -10,13 +10,13 @@
     <CellData>
     </CellData>
     <Coordinates>
-      <DataArray type="Float64" Name="Array 0x1306073f0" format="binary" RangeMin="0" RangeMax="6000">
+      <DataArray type="Float64" Name="Array 0x12d0b2f30" format="binary" RangeMin="0" RangeMax="6000">
         AQAAAACAAAAoAwAAFwEAAA==eJxNzi1MQmEUxvEbCQYiwUAkGAgG4g0GEiMYSI7AHIndOeacwd0x5phzehVUBJUroCL4cf2EQCAaDEQD4UaDgUgwGJ5/8C2/ved9zzmPYfw/CVOuSHuVe06O1qivS3OT9y3p2/wrSLfI/22ZLtG3I8O79O/Rv4+OHB8w71B6ZeZWpHPE/GNpnbCnKpOn7KvJaJ29ZzJ4zv4LOUW/IbMud0xdkgvjTfJhrEVOjLTJi6ErcmPgmvw4Q+tGfmO6I78weSs/0OzKAUZ7soPhO1nF4L0sofEgN3CK2Ufp44JHHdvo4/yTTGEFxzj3LONYxBH+YuxF5tHDH4y8ygy6OMHQm1xGBz8x8C6X0MYhznCxL62++QcopYc0
       </DataArray>
-      <DataArray type="Float64" Name="Array 0x130607530" format="binary" RangeMin="0" RangeMax="0">
+      <DataArray type="Float64" Name="Array 0x12d0b1780" format="binary" RangeMin="0" RangeMax="0">
         AQAAAACAAAAIAAAACwAAAA==eJxjYIAAAAAIAAE=
       </DataArray>
-      <DataArray type="Float64" Name="Array 0x130608990" format="binary" RangeMin="0" RangeMax="3000">
+      <DataArray type="Float64" Name="Array 0x12d09e3c0" format="binary" RangeMin="0" RangeMax="3000">
         AQAAAACAAACYAQAAqQAAAA==eJxNzikOwkAYhuFKJLISiUQiR6JIJYqgCIo0CASCNIQQQggtW9lpS1ksR+AIHKFHQCIR3ysY82Rm/uWzrP9TNbIuvSb3tnx1eO9K0+O/LzOPuoGMh9SPZGNM30QWpvTP6PcxkO858xbyuWTuSgZr5ofS3bBnK50d+/aydGDvUeZP7D/LD2aRbMXcsZaQCysX8mE5JScWr+RF+0ZuzN3Jj190H+YHwvo/dQ==
       </DataArray>
     </Coordinates>

diff --git a/tests/files/sol_analytique_weakly_an_2d_tt.vtr b/tests/files/sol_analytique_weakly_an_2d_tt.vtr