Merge pull request #1203 from danielpeter/devel

.gitattributes

@@ -7,3 +7,4 @@
 *.semc binary
 *.adj binary
 *.ascii binary
+*.bin binary

.github/workflows/CI.yml

@@ -574,4 +574,28 @@ jobs:
       run: ./.github/scripts/run_tests.sh
       shell: bash
 
+  linuxTest_12:
+    name: Test run example 11 - Marmousi2
+    runs-on: ubuntu-latest
+    needs: [linuxCheck]
+
+    steps:
+    - uses: actions/checkout@v3
+
+    - name: Install packages
+      run: ./.github/scripts/run_install.sh
+      shell: bash
+
+    - name: Run build
+      env:
+        TESTFLAGS: --with-mpi --enable-vectorization
+      run: ./.github/scripts/run_build.sh
+      shell: bash
+
+    - name: Run test
+      env:
+        TESTDIR: EXAMPLES/Marmousi2
+      run: ./.github/scripts/run_tests.sh
+      shell: bash
+
 

.travis.yml

@@ -124,6 +124,9 @@ env:
     # infinite_homogeneous_moment_tensor_vertical_dip_slip example
     - TESTID=29 TESTDIR=21 TESTCOV=0 TESTFLAGS="" CUDA=false
 
+    # Marmousi2
+    - TESTID=30 TESTDIR=22 TESTCOV=0 TESTFLAGS="--enable-vectorization --with-mpi" CUDA=false
+
 
 # additional ARM tests
 jobs:

.travis/run_tests.sh

@@ -40,6 +40,7 @@ case "$TESTDIR" in
 19) dir=EXAMPLES/moving_sources_acoustic/ ;;
 20) dir=EXAMPLES/anisotropic_isotropic_model/ ;;
 21) dir=EXAMPLES/infinite_homogeneous_moment_tensor_vertical_dip_slip/ ;;
+22) dir=EXAMPLES/Marmousi2/ ;;
 *) dir=EXAMPLES/simple_topography_and_also_a_simple_fluid_layer/ ;;
 esac
 
@@ -388,9 +389,9 @@ if [ "$TESTCOV" == "1" ] && [ "$TESTID" == "1" ]; then
   ## elastic kernel example Tromp2005_kernel/ w/ NO_BACKWARD_RECONSTRUCTION
   ##
   echo "##################################################################"
-  echo "EXAMPLES/Tromp2005_kernel/"
+  echo "EXAMPLES/Tromp2005_kernel"
   echo
-  cd EXAMPLES/Tromp2005_kernel/
+  cd EXAMPLES/Tromp2005_kernel
   sed -i "s:^NSTEP .*:NSTEP    = 500:" DATA/Par_file
   sed -i "s:^NO_BACKWARD_RECONSTRUCTION .*:NO_BACKWARD_RECONSTRUCTION = .true.:" DATA/Par_file
   sed -i "s:^NTSTEP_BETWEEN_COMPUTE_KERNELS .*:NTSTEP_BETWEEN_COMPUTE_KERNELS = 12:" DATA/Par_file
@@ -402,6 +403,22 @@ if [ "$TESTCOV" == "1" ] && [ "$TESTID" == "1" ]; then
 fi
 echo -en 'travis_fold:end:coverage.no_backward\\r'
 
+echo 'Coverage...' && echo -en 'travis_fold:start:coverage.marmousi2\\r'
+if [ "$TESTCOV" == "1" ] && [ "$TESTID" == "1" ]; then
+  ##
+  ## Marmousi2
+  ##
+  echo "##################################################################"
+  echo "EXAMPLES/Marmousi2"
+  echo
+  cd EXAMPLES/Marmousi2
+  sed -i "s:^NSTEP .*:NSTEP    = 10:" DATA/Par_file
+  ./run_this_example.sh
+  if [[ $? -ne 0 ]]; then exit 1; fi
+  # only for coverage, comparison would fail: my_test
+  cd $WORKDIR
+fi
+echo -en 'travis_fold:end:coverage.marmousi2\\r'
 
 
 #################################################################

EXAMPLES/Marmousi2_to_be_done_one_day_not_ready_yet/DATA/Par_file → EXAMPLES/Marmousi2/DATA/Par_file

@@ -5,7 +5,7 @@
 #-----------------------------------------------------------
 
 # title of job
-title                           = curved interfaces model w/ compaction gradient (marmousi type)
+title                           = Marmousi2 gridded tomo
 
 # forward or adjoint simulation
 # 1 = forward, 2 = adjoint, 3 = both simultaneously
@@ -17,14 +17,14 @@ NOISE_TOMOGRAPHY                = 0
 SAVE_FORWARD                    = .false.
 
 # parameters concerning partitioning
-NPROC                           = 2              # number of processes
+NPROC                           = 4              # number of processes
 
 # time step parameters
 # total number of time steps
-NSTEP                           = 4000
+NSTEP                           = 5000
 
 # duration of a time step (see section "How to choose the time step" of the manual for how to do this)
-DT                              = 2.5d-4
+DT                              = 1.d-3
 
 # time stepping
 # 1 = Newmark (2nd order), 2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta), 3 = classical RK4 4th-order 4-stage Runge-Kutta
@@ -62,15 +62,15 @@ setup_with_binary_database      = 0
 #   external      - define model using define_external_model subroutine
 #   gll           - read GLL model from binary database file
 #   legacy        - read model from model_velocity.dat_input
-MODEL                           = marmousi
+MODEL                           = default
 
 # Output the model with the requested type, does not save if turn to default or .false.
 # (available output formats: ascii,binary,gll,legacy)
-SAVE_MODEL                      = default
+SAVE_MODEL                      = binary
 
 # external tomography file
 # (used for tomography materials with negative material ids and/or MODEL==tomo settings)
-TOMOGRAPHY_FILE                 = ./DATA/tomo_file.xyz
+TOMOGRAPHY_FILE                 = ./DATA/tomography_model_marmousi2.xyz.bin
 
 #-----------------------------------------------------------
 #
@@ -181,11 +181,11 @@ rec_normal_to_surface           = .false.        # base anglerec normal to surfa
 
 # first receiver set (repeat these 6 lines and adjust nreceiversets accordingly)
 nrec                            = 11             # number of receivers
-xdeb                            = 300.           # first receiver x in meters
-zdeb                            = 2200.          # first receiver z in meters
-xfin                            = 3700.          # last receiver x in meters (ignored if only one receiver)
-zfin                            = 2200.          # last receiver z in meters (ignored if only one receiver)
-record_at_surface_same_vertical = .true.         # receivers inside the medium or at the surface (z values are ignored if this is set to true, they are replaced with the topography height)
+xdeb                            = 1000.          # first receiver x in meters
+zdeb                            = 3450.          # first receiver z in meters
+xfin                            = 11000.         # last receiver x in meters (ignored if only one receiver)
+zfin                            = 3450.          # last receiver z in meters (ignored if only one receiver)
+record_at_surface_same_vertical = .false.        # receivers inside the medium or at the surface (z values are ignored if this is set to true, they are replaced with the topography height)
 
 
 #-----------------------------------------------------------
@@ -211,7 +211,7 @@ APPROXIMATE_HESS_KL             = .false.
 
 # Perfectly Matched Layer (PML) boundaries
 # absorbing boundary active or not
-PML_BOUNDARY_CONDITIONS         = .true.
+PML_BOUNDARY_CONDITIONS         = .false.
 NELEM_PML_THICKNESS             = 3
 ROTATE_PML_ACTIVATE             = .false.
 ROTATE_PML_ANGLE                = 30.
@@ -226,7 +226,7 @@ damping_change_factor_elastic   = 1.0d0
 PML_PARAMETER_ADJUSTMENT        = .false.
 
 # Stacey ABC
-STACEY_ABSORBING_CONDITIONS     = .false.
+STACEY_ABSORBING_CONDITIONS     = .true.
 
 # periodic boundaries
 ADD_PERIODIC_CONDITIONS         = .false.
@@ -268,7 +268,7 @@ tangential_detection_curve_file = dummy                                   # file
 
 # material properties
 # number of model materials
-nbmodels                        = 4
+nbmodels                        = 2
 # available material types (see user manual for more information)
 #   acoustic:              model_number 1 rho Vp 0  0 0 QKappa 9999 0 0 0 0 0 0 (for QKappa use 9999 to ignore it)
 #   elastic:               model_number 1 rho Vp Vs 0 0 QKappa Qmu  0 0 0 0 0 0 (for QKappa and Qmu use 9999 to ignore them)
@@ -286,18 +286,16 @@ nbmodels                        = 4
 #       Please also note that Qmu is always equal to Qs, but Qkappa is in general not equal to Qp.
 #       To convert one to the other see doc/Qkappa_Qmu_versus_Qp_Qs_relationship_in_2D_plane_strain.pdf and
 #       utils/attenuation/conversion_from_Qkappa_Qmu_to_Qp_Qs_from_Dahlen_Tromp_959_960.f90.
-1 1 2700.d0 3000.d0 1732.051d0 0.2 0 9999 9999 0 0 0 0 0 0
-2 1 2500.d0 2700.d0 1600.0     0.2 0 9999 9999 0 0 0 0 0 0
-3 1 2200.d0 2500.d0 1443.375d0 0.2 0 9999 9999 0 0 0 0 0 0
-4 1 2200.d0 2200.d0 1343.375d0 0.2 0 9999 9999 0 0 0 0 0 0
+1 -1   0.d0    0.d0    1.d0    0 0    0    0 0 0 0 0 0 0     # solid from tomo file
+2 1 1010.d0 1500.d0    0.d0    0 0 9999 9999 0 0 0 0 0 0     # top water layer
 
 # file containing interfaces for internal mesh
-interfacesfile                  = interfaces_simple_topo_curved.dat
+interfacesfile                  = interfaces.dat
 
 # geometry of the model (origin lower-left corner = 0,0) and mesh description
 xmin                            = 0.d0           # abscissa of left side of the model
-xmax                            = 4000.d0        # abscissa of right side of the model
-nx                              = 400            # number of elements along X
+xmax                            = 17000.d0       # abscissa of right side of the model
+nx                              = 340            # number of elements along X
 
 # absorbing boundary parameters (see absorbing_conditions above)
 absorbbottom                    = .true.
@@ -306,12 +304,10 @@ absorbtop                       = .false.
 absorbleft                      = .true.
 
 # define the different regions of the model in the (nx,nz) spectral-element mesh
-nbregions                       = 4              # then set below the different regions and model number for each region
+nbregions                       = 2              # then set below the different regions and model number for each region
 # format of each line: nxmin nxmax nzmin nzmax material_number
-1 400  1 100 1
-1 400 101 200 2
-1 400 201 300 3
-300 350 101 200 4
+1 340  1 61 1                                    # solid
+1 340 62 70 2                                    # top water layer
 
 #-----------------------------------------------------------
 #
@@ -321,7 +317,7 @@ nbregions                       = 4              # then set below the different
 
 # interval at which we output time step info and max of norm of displacement
 # (every how many time steps we display information about the simulation. costly, do not use a very small value)
-NTSTEP_BETWEEN_OUTPUT_INFO      = 100
+NTSTEP_BETWEEN_OUTPUT_INFO      = 500
 
 # meshing output
 output_grid_Gnuplot             = .false.        # generate a GNUPLOT file containing the grid, and a script to plot it
@@ -346,7 +342,7 @@ COMPUTE_INTEGRATED_ENERGY_FIELD = .false.
 
 # every how many time steps we draw JPEG or PostScript pictures of the simulation
 # and/or we dump results of the simulation as ASCII or binary files (costly, do not use a very small value)
-NTSTEP_BETWEEN_OUTPUT_IMAGES    = 100
+NTSTEP_BETWEEN_OUTPUT_IMAGES    = 500
 
 # minimum amplitude kept in % for the JPEG and PostScript snapshots; amplitudes below that are muted
 cutsnaps                        = 1.
@@ -363,7 +359,7 @@ DRAW_WATER_IN_BLUE              = .true.         # display acoustic layers as co
 USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step (for instance to create movies in an easier way later)
 
 #### for PostScript snapshots ####
-output_postscript_snapshot      = .true.         # output Postscript snapshot of the results every NTSTEP_BETWEEN_OUTPUT_IMAGES time steps or not
+output_postscript_snapshot      = .false.         # output Postscript snapshot of the results every NTSTEP_BETWEEN_OUTPUT_IMAGES time steps or not
 imagetype_postscript            = 1              # display 1=displ vector 2=veloc vector 3=accel vector; small arrows are displayed for the vectors
 meshvect                        = .true.         # display mesh on PostScript plots or not
 modelvect                       = .false.        # display velocity model on PostScript plots or not

EXAMPLES/Marmousi2_to_be_done_one_day_not_ready_yet/DATA/SOURCE → EXAMPLES/Marmousi2/DATA/SOURCE

@@ -1,7 +1,7 @@
 ## Source 1
 source_surf                     = .false.        # source inside the medium, or source automatically moved exactly at the surface by the solver
-xs                              = 2500.          # source location x in meters
-zs                              = 2500.          # source location z in meters (zs is ignored if source_surf is set to true, it is replaced with the topography height)
+xs                              = 5000.          # source location x in meters
+zs                              = 3450.          # source location z in meters (zs is ignored if source_surf is set to true, it is replaced with the topography height)
 ## Source type parameters:
 #  1 = elastic force or acoustic pressure
 #  2 = moment tensor
@@ -33,13 +33,13 @@ source_type                     = 1
 #  9 = burst,
 # 10 = Sinus source time function,
 # 11 = Marmousi Ormsby wavelet
-time_function_type              = 11
+time_function_type              = 1
 # If time_function_type == 8, enter below the custom source file to read (two columns file with time and amplitude) :
 # (For the moment dt must be equal to the dt of the simulation. File name cannot exceed 150 characters)
 # IMPORTANT: do NOT put quote signs around the file name, just put the file name itself otherwise the run will stop
 name_of_source_file             = YYYYYYYYYYYYYYYYYY   # Only for option 8 : file containing the source wavelet
 burst_band_width                = 0.             # Only for option 9 : band width of the burst
-f0                              = 35.0           # dominant source frequency (Hz) if not Dirac or Heaviside
+f0                              = 5.0            # dominant source frequency (Hz) if not Dirac or Heaviside
 tshift                          = 0.0            # time shift when multi sources (if one source, must be zero)
 ## Force source
 # angle of the source (for a force only); for a plane wave, this is the incidence angle; for moment tensor sources this is unused

EXAMPLES/Marmousi2/DATA/interfaces.dat

@@ -0,0 +1,25 @@
+# number of interfaces
+3 
+#
+# for each interface below, we give the number of points and then x,z for each point
+#
+# interface number 1 (bottom of the mesh)
+ 2
+     0.d0     0.d0
+ 17000.d0     0.d0
+# interface number 2 (water bottom)
+ 2
+    0.d0   3050.d0
+ 17000.d0  3050.d0
+
+# interface number 3 (top)
+ 2
+    0.d0   3500.d0
+ 17000.d0  3500.d0
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+# layer number 1 (bottom layer)
+ 61
+# layer number 2 (water layer)
+  9

EXAMPLES/Marmousi2/DATA/tomography_model_marmousi2.xyz.bin.info

@@ -0,0 +1,26 @@
+# tomography model - converted using script interpolate_Marmousi_2_tomo.py
+#
+# providence
+# created             : 2024-02-17 13:05:58.535336
+# command             : ./interpolate_Marmousi_2_tomo.py --smooth=5.0 --replace-water --binary
+#
+# model format
+# model type          : Marmousi2
+#                       Gaussian smoothed for maximum frequency: 5.0 (Hz)
+#                       using solid domain for water layer in Marmousi2
+# coordinate format   : x / z  # z-direction (positive up)
+#
+#origin_x #origin_z #end_x #end_z
+0.0 0.0 17000.0 3500.0
+#dx #dz
+10.0 10.0
+#nx #nz
+1701 351
+#vp_min #vp_max #vs_min #vs_max #density_min #density_max
+1091.5262451171873 4699.970703125 291.0799255371093 2751.93212890625 1630.2421874999998 2626.977294921875
+#
+# data records - format:
+#x #z #vp #vs #density
+#
+# data file: ./DATA/tomography_model_marmousi2.xyz.bin
+#

EXAMPLES/Marmousi2/README.md

@@ -0,0 +1,58 @@
+Readme
+======
+
+This example is demonstrating how to use Marmousi2 based on a tomography file.
+The meshing is done by the internal mesher, using a regular mesh with an element size of ~50m.
+You can change this mesh by modifying the corresponding section in the `DATA/Par_file`.
+
+
+## Tomography file
+
+  We provide a Marmousi2 tomography file in `DATA/tomography_model_marmousi2.xyz.bin` that was created
+  using the script `interpolate_Marmousi_2_tomo.py`:
+  ```
+  ./interpolate_Marmousi_2_tomo.py --smooth=5.0 --replace-water --binary
+  ```
+
+  This script downloads the original, gridded Marmousi2 SEGY-files and together with the above options, it smooths the original model
+  with a Gaussian kernel corresponding to a 5 Hz signal wavelength. Furthermore, it replaces the top water layer with the solid velocities
+  of the upper-most elastic layer. This avoids having problems with reading in the tomography file in SPECFEM2D as well as smearing out
+  the fluid-solid interface by the interpolation of the tomography model velocities onto the mesh.
+
+
+## Mesh setup
+
+  The default mesh setup in the `DATA/Par_file` creates two layers, the top water layer and the solid elastic layer below.
+  For the water layer, the velocities are set to the original water velocities from Marmousi2.
+  The solid, elastic layer below defines a tomographic material. This will interpolate the tomography file velocities from Marmousi2
+  onto all the GLL points in this solid layer.
+
+  In case you want the water layer removed, you can choose:
+  - define the whole mesh as a single solid layer with a tomographic material. This will interpolate the provided tomography file, which has
+    the water replaced by the top-most elastic velocities, onto the whole region.
+    Thus, you would only need to modify in `DATA/Par_file` the `nbmodels` section.
+
+  - adapt the mesh position such that it falls within the solid domain. Marmousi2 has a top water layer with a depth of 450 m.
+    Thus, you can specify the mesh to lie within a vertical range of 0 to 3050 m to be in the solid part only.
+    For this, you would modify in `DATA/Par_file` the `nbmodels` and `nbregions` sections, as well as changing
+    the `DATA/interfaces.dat` file accordingly.
+
+There are many ways how to run external models, and in particular how to interpolate external model velocities onto a SPECFEM2D grid.
+This example demonstrates the use of a tomography file to achieve this, but might not be the best option for your own use cases.
+
+Please consider contributing more examples to this package!
+
+
+## References
+
+Martin, G. S., 2004,
+*The Marmousi2 model, elastic synthetic data, and an analysis of imaging and AVO in a structurally complex environment*:
+Master’s thesis. University of Houston. Retrieved from www.agl.uh.edu/pdf/theses/martin.pdf
+
+Martin, G. S., Marfurt, K. J., and Larsen, S., 2005,
+*Marmousi‐2: An updated model for the investigation of AVO in structurally complex areas*:
+SEG Technical Program Expanded Abstracts 2002, 1979–1982. doi:10.1190/1.1817083
+
+Martin, G. S., Wiley, R., and Marfurt, K. J., 2006,
+*Marmousi2: An elastic upgrade for Marmousi*:
+The Leading Edge, 25, 156–166. doi:10.1190/1.2172306