Merge pull request #1246 from danielpeter/devel

.travis.yml

@@ -17,6 +17,10 @@ compiler: gcc
 # deprecated
 #sudo: required
 
+# turns off submodule fetching by default
+git:
+  submodules: false
+
 env:
   global:
     - FC=gfortran

README.md

@@ -60,6 +60,11 @@ involved in the community and keep them in the specfem3d github wiki:
 [specfem3d wiki](https://github.com/SPECFEM/specfem3d/wiki)
 
 
+## Our contributors :sparkles:
+
+[![SPECFEM2D contributors](https://contrib.rocks/image?repo=SPECFEM/specfem2d)](https://github.com/SPECFEM/specfem2d/graphs/contributors)
+
+
 ## Computational Infrastructure for Geodynamics (CIG)
 
 SPECFEM2D is part of the software that is hosted by the Computational Infrastructure for Geodynamics (CIG). It is available on the CIG website [here (SPECFEM2D)](https://geodynamics.org/resources/specfem2d).

src/meshfem2D/repartition_coupling.f90

@@ -433,6 +433,7 @@ subroutine periodic_edges_repartitioning(elmnts_l,nnodes,nodes_coords,PERIODIC_H
   ! user output
   write(IMAIN,*) 'done detecting points for periodic boundary conditions.'
   write(IMAIN,*) 'number of periodic elements found and grouped in the same partition: ',count(is_periodic)
+  write(IMAIN,*)
   call flush_IMAIN()
 
   ! loop on all the elements to find the first partition that contains a periodic element

src/specfem2D/plot_post.F90

@@ -57,7 +57,8 @@ subroutine plot_post()
                          fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge,num_fluid_solid_edges, &
                          fluid_poro_acoustic_ispec,fluid_poro_acoustic_iedge,num_fluid_poro_edges, &
                          solid_poro_poroelastic_ispec,solid_poro_poroelastic_iedge,num_solid_poro_edges, &
-                         myrank,NPROC
+                         myrank,NPROC, &
+                         P_SV
 
   use shared_parameters, only: subsamp_postscript,imagetype_postscript,interpol, &
     meshvect,modelvect, &
@@ -192,7 +193,12 @@ subroutine plot_post()
   ratio_page = min(rpercentz*sizez/(zmax-zmin),rpercentx*sizex/(xmax-xmin)) / 100.d0
 
   ! compute the maximum of the norm of the vector
-  dispmax = maxval(sqrt(vector_field_display(1,:)**2 + vector_field_display(2,:)**2))
+  if (P_SV) then
+    dispmax = maxval(sqrt(vector_field_display(1,:)**2 + vector_field_display(2,:)**2))
+  else
+    ! SH (membrane) waves, plot y-component
+    dispmax = maxval(abs(vector_field_display(1,:)))
+  endif
 
   call max_all_all_dp(dispmax, dispmax_glob)
   dispmax = dispmax_glob
@@ -1516,22 +1522,39 @@ subroutine plot_post()
           Uxinterp(i,j) = 0.d0
           Uzinterp(i,j) = 0.d0
 
-          do k = 1,NGLLX
-            do l= 1,NGLLX
-              if (AXISYM) then
-                if (is_on_the_axis(ispec)) then
-                  Uxinterp(i,j) = Uxinterp(i,j) + vector_field_display(1,ibool(k,l,ispec))*flagrange_GLJ(k,i)*flagrange_GLJ(l,j)
-                  Uzinterp(i,j) = Uzinterp(i,j) + vector_field_display(2,ibool(k,l,ispec))*flagrange_GLJ(k,i)*flagrange_GLJ(l,j)
+          if (P_SV) then
+            do k = 1,NGLLX
+              do l= 1,NGLLX
+                if (AXISYM) then
+                  if (is_on_the_axis(ispec)) then
+                    Uxinterp(i,j) = Uxinterp(i,j) + vector_field_display(1,ibool(k,l,ispec))*flagrange_GLJ(k,i)*flagrange_GLJ(l,j)
+                    Uzinterp(i,j) = Uzinterp(i,j) + vector_field_display(2,ibool(k,l,ispec))*flagrange_GLJ(k,i)*flagrange_GLJ(l,j)
+                  else
+                    Uxinterp(i,j) = Uxinterp(i,j) + vector_field_display(1,ibool(k,l,ispec))*flagrange(k,i)*flagrange(l,j)
+                    Uzinterp(i,j) = Uzinterp(i,j) + vector_field_display(2,ibool(k,l,ispec))*flagrange(k,i)*flagrange(l,j)
+                  endif
                 else
                   Uxinterp(i,j) = Uxinterp(i,j) + vector_field_display(1,ibool(k,l,ispec))*flagrange(k,i)*flagrange(l,j)
                   Uzinterp(i,j) = Uzinterp(i,j) + vector_field_display(2,ibool(k,l,ispec))*flagrange(k,i)*flagrange(l,j)
                 endif
-              else
-                Uxinterp(i,j) = Uxinterp(i,j) + vector_field_display(1,ibool(k,l,ispec))*flagrange(k,i)*flagrange(l,j)
-                Uzinterp(i,j) = Uzinterp(i,j) + vector_field_display(2,ibool(k,l,ispec))*flagrange(k,i)*flagrange(l,j)
-              endif
+              enddo
             enddo
-          enddo
+          else
+            ! SH (membrane) waves, plot y-component
+            do k = 1,NGLLX
+              do l= 1,NGLLX
+                if (AXISYM) then
+                  if (is_on_the_axis(ispec)) then
+                    Uxinterp(i,j) = Uxinterp(i,j) + vector_field_display(1,ibool(k,l,ispec))*flagrange_GLJ(k,i)*flagrange_GLJ(l,j)
+                  else
+                    Uxinterp(i,j) = Uxinterp(i,j) + vector_field_display(1,ibool(k,l,ispec))*flagrange(k,i)*flagrange(l,j)
+                  endif
+                else
+                  Uxinterp(i,j) = Uxinterp(i,j) + vector_field_display(1,ibool(k,l,ispec))*flagrange(k,i)*flagrange(l,j)
+                endif
+              enddo
+            enddo
+          endif
 
           x1 =(xinterp(i,j)-xmin)*ratio_page
           z1 =(zinterp(i,j)-zmin)*ratio_page
@@ -1672,8 +1695,14 @@ subroutine plot_post()
       z1 =(coord(2,ipoin)-zmin)*ratio_page
 
       if (dispmax > 0.d0) then
-        x2 = vector_field_display(1,ipoin)*sizemax_arrows/dispmax
-        z2 = vector_field_display(2,ipoin)*sizemax_arrows/dispmax
+        if (P_SV) then
+          x2 = vector_field_display(1,ipoin)*sizemax_arrows/dispmax
+          z2 = vector_field_display(2,ipoin)*sizemax_arrows/dispmax
+        else
+          ! SH (membrane) waves, plot y-component
+          x2 = vector_field_display(1,ipoin)*sizemax_arrows/dispmax
+          z2 = 0.d0
+        endif
       else
         x2 = 0.d0
         z2 = 0.d0

src/specfem2D/setup_mesh.F90

@@ -134,6 +134,7 @@ subroutine setup_mesh_numbering()
 #else
     write(IMAIN,*) 'Exact total number of grid points in the mesh: ',nglob_total
 #endif
+    write(IMAIN,*)
 
     ! percentage of elements with 2 degrees of freedom per point
     ratio_2DOFs = (nspec_total - nspec_acoustic_total) / dble(nspec_total)
@@ -145,7 +146,6 @@ subroutine setup_mesh_numbering()
     nb_elastic_DOFs  = nint(nglob_total*ratio_2DOFs*2)
 
     if (P_SV) then
-      write(IMAIN,*)
       write(IMAIN,*) 'Approximate number of acoustic degrees of freedom in the mesh: ',nb_acoustic_DOFs
       write(IMAIN,*) 'Approximate number of elastic degrees of freedom in the mesh: ',nb_elastic_DOFs
       write(IMAIN,*) '  (there are 2 degrees of freedom per point for elastic elements)'
@@ -350,7 +350,7 @@ subroutine setup_mesh_periodic_edges()
   ! local parameters
   integer :: ispec,i,j,iglob,iglob2,ier
   double precision :: xmaxval,xminval,ymaxval,yminval,xtol,xtypdist
-  integer :: counter
+  integer :: counter,counter_all
 
 ! allocate an array to make sure that an acoustic free surface is not enforced on periodic edges
   allocate(this_ibool_is_a_periodic_edge(NGLOB),stat=ier)
@@ -362,19 +362,20 @@ subroutine setup_mesh_periodic_edges()
   if (ADD_PERIODIC_CONDITIONS) then
     ! user output
     if (myrank == 0) then
-      write(IMAIN,*)
-      write(IMAIN,*) 'implementing periodic boundary conditions'
-      write(IMAIN,*) 'in the horizontal direction with a periodicity distance of ',PERIODIC_HORIZ_DIST,' m'
+      write(IMAIN,*) 'Periodic boundary conditions:'
+      write(IMAIN,*) '  implementing periodic boundary conditions'
+      write(IMAIN,*) '  in the horizontal direction with a periodicity distance of ',sngl(PERIODIC_HORIZ_DIST),' m'
       if (PERIODIC_HORIZ_DIST <= 0.d0) call stop_the_code( &
 'PERIODIC_HORIZ_DIST should be greater than zero when using ADD_PERIODIC_CONDITIONS')
       write(IMAIN,*)
-      write(IMAIN,*) '*****************************************************************'
-      write(IMAIN,*) '*****************************************************************'
-      write(IMAIN,*) '**** BEWARE: because of periodic conditions, values computed ****'
-      write(IMAIN,*) '****         by check_grid() below will not be reliable       ****'
-      write(IMAIN,*) '*****************************************************************'
-      write(IMAIN,*) '*****************************************************************'
+      write(IMAIN,*) '  *****************************************************************'
+      write(IMAIN,*) '  *****************************************************************'
+      write(IMAIN,*) '  **** BEWARE: because of periodic conditions, values computed ****'
+      write(IMAIN,*) '  ****         by check_grid() below will not be reliable      ****'
+      write(IMAIN,*) '  *****************************************************************'
+      write(IMAIN,*) '  *****************************************************************'
       write(IMAIN,*)
+      call flush_IMAIN()
     endif
 
     ! set up a local geometric tolerance
@@ -419,8 +420,9 @@ subroutine setup_mesh_periodic_edges()
 ! (as implemented in routine createnum_fast() elsewhere in the code). This could be done one day if needed instead
 ! of the very simple double loop below.
     if (myrank == 0) then
-      write(IMAIN,*) 'start detecting points for periodic boundary conditions '// &
+      write(IMAIN,*) '  start detecting points for periodic boundary conditions '// &
                      '(the current algorithm can be slow and could be improved)...'
+      call flush_IMAIN()
     endif
 
     counter = 0
@@ -447,9 +449,31 @@ subroutine setup_mesh_periodic_edges()
       enddo
     enddo
 
-    if (myrank == 0) write(IMAIN,*) 'done detecting points for periodic boundary conditions.'
+    if (myrank == 0) then
+      write(IMAIN,*) '  done detecting points for periodic boundary conditions.'
+      write(IMAIN,*)
+      call flush_IMAIN()
+    endif
 
-    if (counter > 0) write(IMAIN,*) 'implemented periodic conditions on ',counter,' grid points on proc ',myrank
+    if (counter > 0) then
+      write(IMAIN,*) '  implemented periodic conditions on ',counter,' grid points on proc ',myrank
+    endif
+
+    ! check if any points found
+    call sum_all_i(counter,counter_all)
+
+    if (myrank == 0) then
+      write(IMAIN,*) '  total number of grid points found for periodic conditions = ',counter_all
+      write(IMAIN,*)
+      if (counter_all == 0) then
+        write(IMAIN,*) '  No grid points found for periodic conditions.'
+        write(IMAIN,*) '  Detection uses a typical element size ',xtypdist,'and position tolerance ',xtol
+        write(IMAIN,*) '  Please check if periodic horizontal distance ',sngl(PERIODIC_HORIZ_DIST), &
+                       'is coherent with mesh dimensions'
+        write(IMAIN,*)
+      endif
+      call flush_IMAIN()
+    endif
 
   endif ! of if (ADD_PERIODIC_CONDITIONS)
 

src/specfem2D/write_postscript_snapshot.f90

@@ -52,31 +52,31 @@ subroutine write_postscript_snapshot()
   endif
 
   ! determines postscript output type
-  if (P_SV) then
-    select case (imagetype_postscript)
-    case (1)
-      ! displacement
-      if (myrank == 0) write(IMAIN,*) 'drawing displacement vector as small arrows...'
-      call compute_vector_whole_medium(potential_acoustic,displ_elastic,displs_poroelastic)
-    case (2)
-      ! velocity
-      if (myrank == 0) write(IMAIN,*) 'drawing velocity vector as small arrows...'
-      call compute_vector_whole_medium(potential_dot_acoustic,veloc_elastic,velocs_poroelastic)
-    case (3)
-      ! acceleration
-      if (myrank == 0) write(IMAIN,*) 'drawing acceleration vector as small arrows...'
-      call compute_vector_whole_medium(potential_dot_dot_acoustic,accel_elastic,accels_poroelastic)
-    case default
-      call exit_MPI(myrank,'wrong type for PostScript snapshots')
-    end select
+  select case (imagetype_postscript)
+  case (1)
+    ! displacement
+    if (myrank == 0) write(IMAIN,*) 'drawing displacement vector as small arrows...'
+    call compute_vector_whole_medium(potential_acoustic,displ_elastic,displs_poroelastic)
+  case (2)
+    ! velocity
+    if (myrank == 0) write(IMAIN,*) 'drawing velocity vector as small arrows...'
+    call compute_vector_whole_medium(potential_dot_acoustic,veloc_elastic,velocs_poroelastic)
+  case (3)
+    ! acceleration
+    if (myrank == 0) write(IMAIN,*) 'drawing acceleration vector as small arrows...'
+    call compute_vector_whole_medium(potential_dot_dot_acoustic,accel_elastic,accels_poroelastic)
+  case default
+    call exit_MPI(myrank,'wrong type for PostScript snapshots')
+  end select
 
-    ! postscript plotting
-    call plot_post()
-
-  else
-    call exit_MPI(myrank,'cannot draw a SH scalar field as a vector plot, turn PostScript plots off')
+  ! info for SH simulations
+  if (.not. P_SV) then
+    write(IMAIN,*) 'drawing a SH scalar field as a vector plot oriented along x-direction'
   endif
 
+  ! postscript plotting
+  call plot_post()
+
   ! user output
   if (myrank == 0) then
     write(IMAIN,*) 'PostScript file written'