need vp

slimgroup · Sep 30, 2024 · ad860c2 · ad860c2
1 parent af00642
commit ad860c2
Show file tree

Hide file tree

Showing 8 changed files with 89 additions and 161 deletions.
diff --git a/.github/workflows/ci-judi.yml b/.github/workflows/ci-judi.yml
@@ -19,7 +19,7 @@ jobs:
     name: JUDI base on Julia ${{ matrix.version }} - ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
     env:
-      DEVITO_ARCH: gcc-11
+      DEVITO_ARCH: gcc-12
       DEVITO_LANGUAGE: "openmp"
       OMP_NUM_THREADS: 4
       GROUP: "JUDI"

diff --git a/.github/workflows/ci-op.yml b/.github/workflows/ci-op.yml
@@ -19,7 +19,7 @@ jobs:
     name: ${{ matrix.op }} on Julia ${{ matrix.version }} - ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
     env:
-      DEVITO_ARCH: gcc-11
+      DEVITO_ARCH: ${{ matrix.cc }}
       DEVITO_LANGUAGE: "openmp"
       DEVITO_LOGGING: "INFO"
       OMP_NUM_THREADS: ${{ matrix.omp }}
@@ -33,32 +33,38 @@ jobs:
         op: ["ISO_OP", "ISO_OP_FS", "TTI_OP", "TTI_OP_FS", "BASICS"]
         version: ['1']
         omp: [2]
+        cc: ['gcc-11']
 
         include:
           - os: macos-13
             version: '1.6'
             op: "ISO_OP"
             omp: 1
+            cc: gcc-13
 
           - os: macos-13
             version: '1.8'
             op: "ISO_OP"
             omp: 1
+            cc: gcc-13
 
           - os: macos-13
             version: '1.9'
             op: "ISO_OP"
             omp: 1
+            cc: gcc-13
 
           - os: ubuntu-latest
             version: '1.9'
             op: "VISCO_AC_OP"
             omp: 2
+            cc: gcc-11
 
           - os: ubuntu-latest
             version: '1.10'
             op: "ISO_OP"
             omp: 2
+            cc: gcc-11
 
     steps:
       - name: Checkout JUDI

diff --git a/Project.toml b/Project.toml
@@ -12,6 +12,7 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 JOLI = "bb331ad6-a1cf-11e9-23da-9bcb53c69f6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"

diff --git a/examples/scripts/modeling_basic_elastic.jl b/examples/scripts/modeling_basic_elastic.jl
@@ -29,7 +29,7 @@ nxrec = 120
 nyrec = 100
 xrec = range(50f0, stop=1150f0, length=nxrec)
 yrec = 0f0
-zrec = range(0f0, stop=0f0, length=nxrec)
+zrec = range(10f0, stop=10f0, length=nxrec)
 
 # receiver sampling and recording time
 timeR = 1500f0   # receiver recording time [ms]
@@ -41,7 +41,7 @@ recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)
 # Set up source geometry (cell array with source locations for each shot)
 xsrc = 600f0
 ysrc = 0f0
-zsrc = 0f0
+zsrc = 10f0
 
 # source sampling and number of time steps
 timeS = 1500f0   # source length in [ms]
@@ -57,7 +57,7 @@ wavelet = ricker_wavelet(timeS, dtS, f0)
 ###################################################################################################
 
 # Setup operators
-F = judiModeling(model, srcGeometry, recGeometry)
+F = judiModeling(model, srcGeometry, recGeometry; options=Options(space_order=8, free_surface=true))
 q = judiVector(srcGeometry, wavelet)
 
 # Nonlinear modeling

diff --git a/src/pysource/geom_utils.py b/src/pysource/geom_utils.py
@@ -7,7 +7,7 @@
 try:
     from recipes.utils import mirror_source
 except ImportError:
-    def mirror_source(src):
+    def mirror_source(model, src):
         return src
 
 
@@ -20,12 +20,12 @@ def src_rec(model, u, src_coords=None, rec_coords=None, wavelet=None, nt=None):
             src = wavelet
         else:
             src = PointSource(name="src%s" % namef, grid=model.grid, ntime=nt,
-                              coordinates=src_coords, interpolation='sinc', r=3)
+                              coordinates=src_coords)
             src.data[:] = wavelet.view(np.ndarray) if wavelet is not None else 0.
     rcv = None
     if rec_coords is not None:
         rcv = Receiver(name="rcv%s" % namef, grid=model.grid, ntime=nt,
-                       coordinates=rec_coords, interpolation='sinc', r=3)
+                       coordinates=rec_coords)
     return src, rcv
 
 
@@ -62,28 +62,35 @@ def geom_expr(model, u, src_coords=None, rec_coords=None, wavelet=None, fw=True,
     if not model.is_elastic:
         m = model.m * irho
     dt = model.grid.time_dim.spacing
-    geom_expr = []
     src, rcv = src_rec(model, u, src_coords, rec_coords, wavelet, nt)
     model.__init_abox__(src, rcv, fw=fw)
+
+    geom_expr = []
+    # Source
     if src is not None:
         # Elastic inject into diagonal of stress
         if model.is_elastic:
-            for ud in as_tuple(u)[1].diagonal():
-                geom_expr += src.inject(field=ud.forward, expr=src*dt/irho)
+            c = 1 / model.grid.dim
+            src_eq = src.inject(field=as_tuple(u)[1].forward.diagonal(),
+                                expr=c*src*dt/irho)
+            if model.fs:
+                # Free surface
+                src_eq = mirror_source(model, src_eq)
         else:
             # Acoustic inject into pressure
             u_n = as_tuple(u)[0].forward if fw else as_tuple(u)[0].backward
             src_eq = src.inject(field=u_n, expr=src*dt**2/m)
             if model.fs:
                 # Free surface
                 src_eq = mirror_source(model, src_eq)
-            geom_expr += src_eq
+
+        geom_expr += [src_eq]
     # Setup adjoint wavefield sampling at source locations
     if rcv is not None:
         if model.is_elastic:
-            rec_expr = u[1].trace()
+            geom_expr = u[1].trace() / model.grid.dim
         else:
             rec_expr = u[0] if model.is_tti else u
-        adj_rcv = rcv.interpolate(expr=rec_expr)
-        geom_expr += adj_rcv
+        geom_expr += rcv.interpolate(expr=rec_expr)
+
     return geom_expr