updats compute_arrays_source** routines for moment-tensor and force solutions; updates source location calls

Author: danielpeter

src/shared/read_source_file.f90

@@ -308,7 +308,7 @@ subroutine read_source_file(NSOURCES,BROADCAST_AFTER_READ)
       case (0)
         write(IMAIN,*) '  Normalized Gaussian:'
         write(IMAIN,*) '  Frequency = ',f0_source(i_source), &
-                       ' (-> half-duration = ',sngl(1.d0/f0_source(i_source)),' s)'
+                       ' (half-duration = ',sngl(1.d0/f0_source(i_source)),' s)'
         if (f0_source(i_source) == f0_sampling) write(IMAIN,*) '  (Frequency limited by sampling rate DT ',sngl(DT),' s)'
         write(IMAIN,*) '  delay     = ',tshift_src(i_source)
       case (1)
@@ -349,19 +349,19 @@ subroutine read_source_file(NSOURCES,BROADCAST_AFTER_READ)
       case (12)
         write(IMAIN,*) '  Brune source time function'
         write(IMAIN,*) '  Frequency = ',f0_source(i_source), &
-                       ' (-> rise time = ',sngl(1.d0/f0_source(i_source)),' s)'
+                       ' (rise time = ',sngl(1.d0/f0_source(i_source)),' s)'
         write(IMAIN,*) '  delay     = ',tshift_src(i_source)
       case (13)
         write(IMAIN,*) '  Smoothed Brune source time function'
         write(IMAIN,*) '  Frequency = ',f0_source(i_source), &
-                       ' (-> rise time = ',sngl(1.d0/f0_source(i_source)),' s)'
+                       ' (rise time = ',sngl(1.d0/f0_source(i_source)),' s)'
         write(IMAIN,*) '  delay     = ',tshift_src(i_source)
       case (14)
         write(IMAIN,*) '  Regularized Yoffe:'
-        write(IMAIN,*) '  Frequency        = ',f0_source(i_source),&
-                       ' (-> slip acceleration duration: T_acc = ',sngl(1.d0/f0_source(i_source)),' s)'
+        write(IMAIN,*) '  Frequency        = ',f0_source(i_source), &
+                       ' (slip acceleration duration: T_acc = ',sngl(1.d0/f0_source(i_source)),' s)'
         write(IMAIN,*) '  burst band width = ',burst_band_width(i_source), &
-                       ' (-> effective final duration  : T_eff = ',sngl(1.d0/burst_band_width(i_source)),' s)'
+                       ' (effective final duration  : T_eff = ',sngl(1.d0/burst_band_width(i_source)),' s)'
         write(IMAIN,*) '  delay            = ',tshift_src(i_source)
       case default
         call stop_the_code('Error invalid source time function type! must be between 1 and 9, exiting...')

src/specfem2D/comp_source_time_function.f90

@@ -782,7 +782,7 @@ double precision function comp_source_time_function_Yoffe(t,f0,burst_band_width)
 
   ! regularized Yoffe function defined in Appendix (A13) - (A20) of
   ! Tinti et al. (2005),
-  ! A Kinematic Source-Time Function Compatible with Earthquake Dynamics,
+  ! A kinematic source-time function compatible with earthquake dynamics,
   ! BSSA, 95 (4), 1211-1223. https://doi.org/10.1785/0120040177
 
   ! fixed parameter example:

src/specfem2D/compute_add_sources_acoustic.f90

@@ -98,13 +98,12 @@ subroutine compute_add_sources_acoustic_moving_sources(potential_dot_dot_acousti
 
   use specfem_par, only: ispec_is_acoustic,nglob_acoustic, &
                          NSOURCES,source_type,source_time_function, &
-                         islice_selected_source,ispec_selected_source, &
-                         hxis_store,hgammas_store,ibool,kappastore,myrank,DT,t0,tshift_src, &
+                         islice_selected_source,ispec_selected_source,sourcearrays, &
+                         ibool,kappastore,myrank,DT,t0,tshift_src, &
                          coord,nspec,nglob,xigll,zigll,NPROC,xi_source, &
                          gamma_source,coorg,knods,NGNOD,npgeo,iglob_source,x_source,z_source, &
-                         vx_source,vz_source, time_stepping_scheme, &
-                         SOURCE_IS_MOVING, &
-                         hxis,hpxis,hgammas,hpgammas
+                         anglesource,vx_source,vz_source,time_stepping_scheme, &
+                         AXISYM,xiglj,is_on_the_axis,SOURCE_IS_MOVING
 
   use moving_sources_par, only: locate_source_moving
 
@@ -113,10 +112,14 @@ subroutine compute_add_sources_acoustic_moving_sources(potential_dot_dot_acousti
   real(kind=CUSTOM_REAL), dimension(nglob_acoustic),intent(inout) :: potential_dot_dot_acoustic
   integer,intent(in) :: it,i_stage
 
-  !local variables
+  ! local variables
   integer :: i_source,i,j,iglob,ispec
-  double precision :: hlagrange
+  real(kind=CUSTOM_REAL) :: stf_used
   double precision :: xsrc,zsrc,timeval,t_used
+  ! Lagrange interpolators at source position
+  double precision, dimension(NGLLX) :: hxis,hpxis
+  double precision, dimension(NGLLZ) :: hgammas,hpgammas
+  ! single source array
   real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLZ) :: sourcearray
 
   ! checks if anything to do
@@ -157,6 +160,7 @@ subroutine compute_add_sources_acoustic_moving_sources(potential_dot_dot_acousti
     enddo
   endif
 
+  ! updates source positions and re-calculates source arrays
   do i_source = 1,NSOURCES
     if (abs(source_time_function(i_source,it,i_stage)) > TINYVAL) then
       t_used = (timeval-t0-tshift_src(i_source))
@@ -165,40 +169,55 @@ subroutine compute_add_sources_acoustic_moving_sources(potential_dot_dot_acousti
       xsrc = x_source(i_source) + vx_source(i_source)*t_used
       zsrc = z_source(i_source) + vz_source(i_source)*t_used
 
-      ! collocated force source
+      ! gets source positioning
       ! TODO: this would be more efficient compled with first guess as in init_moving_sources_GPU
       !call locate_source_moving(xsrc,zsrc, &
       !                   ispec_selected_source(i_source),islice_selected_source(i_source), &
-      !                   NPROC,myrank,xi_source(i_source),gamma_source(i_source),.true.)
+      !                   NPROC,myrank,xi_source(i_source),gamma_source(i_source),..
       call locate_source(ibool,coord,nspec,nglob,xigll,zigll, &
                          xsrc,zsrc, &
                          ispec_selected_source(i_source),islice_selected_source(i_source), &
                          NPROC,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,NGNOD,npgeo, &
-                         iglob_source(i_source),.true.)
+                         iglob_source(i_source),source_type(i_source))
 
       ! print *,ispec_selected_source(i_source) > nspec, "xmin:", &
       !               coord(1,ibool(1,1,ispec_selected_source(i_source))), &
       !               "xmax:", coord(1,ibool(NGLLX,1,ispec_selected_source(i_source)))
 
-      ! define and store Lagrange interpolators (hxis,hpxis,hgammas,hpgammas) at all the sources
-      !if (AXISYM) then
-      !  if (is_on_the_axis(ispec_selected_source(i_source)) .and. myrank == islice_selected_source(i_source)) then
-      !    call lagrange_any(xi_source(i_source),NGLJ,xiglj,hxis,hpxis)
-      !    !do j = 1,NGLJ ! ABAB same result with that loop, this is good
-      !    !  hxis(j) = hglj(j-1,xi_source(i),xiglj,NGLJ)
-      !    !enddo
-      !  else
-      !    call lagrange_any(xi_source(i_source),NGLLX,xigll,hxis,hpxis)
-      !  endif
-      !else
-        call lagrange_any(xi_source(i_source),NGLLX,xigll,hxis,hpxis)
-      !endif
-      call lagrange_any(gamma_source(i_source),NGLLZ,zigll,hgammas,hpgammas)
-
-      ! stores Lagrangians for source
-      hxis_store(i_source,:) = hxis(:)
-      hgammas_store(i_source,:) = hgammas(:)
-
+      if (myrank == islice_selected_source(i_source)) then
+        ! element containing source
+        ispec = ispec_selected_source(i_source)
+
+        ! only for acoustic source elements
+        if (.not. ispec_is_acoustic(ispec)) cycle
+
+        ! Lagrange interpolators
+        if (AXISYM) then
+          if (is_on_the_axis(ispec)) then
+            call lagrange_any(xi_source(i_source),NGLJ,xiglj,hxis,hpxis)
+          else
+            call lagrange_any(xi_source(i_source),NGLLX,xigll,hxis,hpxis)
+          endif
+        else
+          call lagrange_any(xi_source(i_source),NGLLX,xigll,hxis,hpxis)
+        endif
+        call lagrange_any(gamma_source(i_source),NGLLZ,zigll,hgammas,hpgammas)
+
+        ! computes source arrays
+        sourcearray(:,:,:) = 0._CUSTOM_REAL
+
+        select case (source_type(i_source))
+        case (1)
+          ! collocated force source
+          call compute_arrays_source_forcesolution(ispec,hxis,hgammas,sourcearray,anglesource(i_source))
+        case (2)
+          ! moment-tensor source
+          call exit_MPI(myrank,'Cannot have moment tensor source in acoustic element')
+        end select
+
+        ! stores sourcearray for all sources
+        sourcearrays(:,:,:,i_source) = sourcearray(:,:,:)
+      endif
     endif
   enddo
 
@@ -212,6 +231,10 @@ subroutine compute_add_sources_acoustic_moving_sources(potential_dot_dot_acousti
       ispec = ispec_selected_source(i_source)
 
       if (ispec_is_acoustic(ispec)) then
+
+        ! source time function
+        stf_used = source_time_function(i_source,it,i_stage)
+
         ! collocated force
         ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
         ! the sign is negative because pressure p = - Chi_dot_dot therefore we need
@@ -222,16 +245,10 @@ subroutine compute_add_sources_acoustic_moving_sources(potential_dot_dot_acousti
             do i = 1,NGLLX
               iglob = ibool(i,j,ispec)
 
-              hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
-              sourcearray(1,i,j) = hlagrange
-
               potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
-                                      + source_time_function(i_source,it,i_stage) * sourcearray(1,i,j) / kappastore(i,j,ispec)
+                                                + sourcearrays(1,i,j,i_source) * stf_used / kappastore(i,j,ispec)
             enddo
           enddo
-          ! moment tensor
-          else if (source_type(i_source) == 2) then
-            call exit_MPI(myrank,'Cannot have moment tensor source in acoustic element')
         endif
       endif
     endif ! if this processor core carries the source and the source element is acoustic

src/specfem2D/compute_add_sources_viscoelastic.f90

@@ -109,9 +109,9 @@ subroutine compute_add_sources_viscoelastic_moving_sources(accel_elastic,it,i_st
                          ibool,coord,nspec,nglob,xigll,zigll,NPROC, &
                          xi_source,gamma_source,coorg,knods,NGNOD,npgeo,iglob_source,x_source,z_source, &
                          vx_source,vz_source,DT,t0,myrank, &
-                         time_stepping_scheme,hxis_store,hgammas_store,tshift_src,source_type,ispec_is_acoustic, &
-                         hxis,hpxis,hgammas,hpgammas,anglesource,ispec_is_poroelastic,Mxx,Mxz,Mzz,gammax,gammaz,xix,xiz, &
-                         AXISYM,xiglj,is_on_the_axis,initialfield,SOURCE_IS_MOVING
+                         time_stepping_scheme,tshift_src,source_type, &
+                         anglesource,Mxx,Mxz,Mzz,gammax,gammaz,xix,xiz, &
+                         AXISYM,xiglj,is_on_the_axis,SOURCE_IS_MOVING
 
   use moving_sources_par, only: locate_source_moving
 
@@ -120,11 +120,13 @@ subroutine compute_add_sources_viscoelastic_moving_sources(accel_elastic,it,i_st
   real(kind=CUSTOM_REAL), dimension(NDIM,nglob_elastic) :: accel_elastic
   integer :: it, i_stage
 
-  !local variable
+  ! local variables
   integer :: i_source,i,j,iglob,ispec
   real(kind=CUSTOM_REAL) :: stf_used
-  double precision :: hlagrange
   double precision :: xsrc,zsrc,timeval,t_used
+  ! Lagrange interpolators at source position
+  double precision, dimension(NGLLX) :: hxis,hpxis
+  double precision, dimension(NGLLZ) :: hgammas,hpgammas
   ! single source array
   real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLZ) :: sourcearray
 
@@ -165,6 +167,7 @@ subroutine compute_add_sources_viscoelastic_moving_sources(accel_elastic,it,i_st
     enddo
   endif
 
+  ! updates source positions and re-calculates source arrays
   do i_source = 1,NSOURCES
     if (abs(source_time_function(i_source,it,i_stage)) > TINYVAL) then
       t_used = (timeval-t0-tshift_src(i_source))
@@ -173,43 +176,28 @@ subroutine compute_add_sources_viscoelastic_moving_sources(accel_elastic,it,i_st
       xsrc = x_source(i_source) + vx_source(i_source)*t_used
       zsrc = z_source(i_source) + vz_source(i_source)*t_used
 
-      ! collocated force source
-      if (source_type(i_source) == 1) then
-        ! TODO: this would be more efficient compled with first guess as in init_moving_sources_GPU()
-        !call locate_source_moving(xsrc,zsrc, &
-        !                   ispec_selected_source(i_source),islice_selected_source(i_source), &
-        !                   NPROC,myrank,xi_source(i_source),gamma_source(i_source),.true.)
-        call locate_source(ibool,coord,nspec,nglob,xigll,zigll, &
-                           xsrc,zsrc, &
-                           ispec_selected_source(i_source),islice_selected_source(i_source), &
-                           NPROC,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,NGNOD,npgeo, &
-                           iglob_source(i_source),.true.)
-
-      else if (source_type(i_source) == 2) then
-        ! moment-tensor source
-        call locate_source(ibool,coord,nspec,nglob,xigll,zigll, &
-                           xsrc,zsrc, &
-                           ispec_selected_source(i_source),islice_selected_source(i_source), &
-                           NPROC,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,NGNOD,npgeo, &
-                           iglob_source(i_source),.false.)
-
-      else if (.not. initialfield) then
-
-        call exit_MPI(myrank,'incorrect source type')
+      ! gets source positioning
+      ! TODO: this would be more efficient compled with first guess as in init_moving_sources_GPU()
+      !call locate_source_moving(xsrc,zsrc, &
+      !                   ispec_selected_source(i_source),islice_selected_source(i_source), &
+      !                   NPROC,myrank,xi_source(i_source),gamma_source(i_source),..
+      call locate_source(ibool,coord,nspec,nglob,xigll,zigll, &
+                         xsrc,zsrc, &
+                         ispec_selected_source(i_source),islice_selected_source(i_source), &
+                         NPROC,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,NGNOD,npgeo, &
+                         iglob_source(i_source),source_type(i_source))
 
-      endif
+      if (myrank == islice_selected_source(i_source)) then
+        ! element containing source
+        ispec = ispec_selected_source(i_source)
 
-      ispec = ispec_selected_source(i_source)
+        ! only for elastic source elements
+        if (.not. ispec_is_elastic(ispec)) cycle
 
-      ! source element is elastic
-      if (ispec_is_elastic(ispec)) then
         ! Lagrange interpolators
         if (AXISYM) then
           if (is_on_the_axis(ispec)) then
             call lagrange_any(xi_source(i_source),NGLJ,xiglj,hxis,hpxis)
-            !do j = 1,NGLJ ! ABAB same result with that loop, this is good
-            !  hxis(j) = hglj(j-1,xi_source(i_source),xiglj,NGLJ)
-            !enddo
           else
             call lagrange_any(xi_source(i_source),NGLLX,xigll,hxis,hpxis)
           endif
@@ -231,79 +219,21 @@ subroutine compute_add_sources_viscoelastic_moving_sources(accel_elastic,it,i_st
           endif
         endif
 
-        ! stores Lagrangians for source
-        hxis_store(i_source,:) = hxis(:)
-        hgammas_store(i_source,:) = hgammas(:)
-
+        ! computes source arrays
         sourcearray(:,:,:) = 0._CUSTOM_REAL
 
-        ! computes source arrays
         select case (source_type(i_source))
         case (1)
           ! collocated force source
-          do j = 1,NGLLZ
-            do i = 1,NGLLX
-              hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
-
-              ! source element is acoustic
-              if (ispec_is_acoustic(ispec)) then
-                sourcearray(:,i,j) = real(hlagrange,kind=CUSTOM_REAL)
-              endif
-
-              ! source element is elastic
-              if (ispec_is_elastic(ispec)) then
-                if (P_SV) then
-                  ! P_SV case
-!                  sourcearray(1,i,j) = real(- sin(anglesource(i_source)) * hlagrange,kind=CUSTOM_REAL)
-!                  sourcearray(2,i,j) = real(  cos(anglesource(i_source)) * hlagrange,kind=CUSTOM_REAL)
-!! DK DK May 2018: the sign of the source was inverted compared to the analytical solution for a simple elastic benchmark
-!! DK DK May 2018: with a force source (the example that is in EXAMPLES/check_absolute_amplitude_of_force_source_seismograms),
-!! DK DK May 2018: which means that the sign was not right here. I changed it. Please do NOT revert that change,
-!! DK DK May 2018: otherwise the code will give inverted seismograms compared to analytical solutions for benchmarks,
-!! DK DK May 2018: and more generally compared to reality
-                  sourcearray(1,i,j) = real(+ sin(anglesource(i_source)) * hlagrange,kind=CUSTOM_REAL)
-                  sourcearray(2,i,j) = real(- cos(anglesource(i_source)) * hlagrange,kind=CUSTOM_REAL)
-                else
-                  ! SH case (membrane)
-                  sourcearray(:,i,j) = real(hlagrange,kind=CUSTOM_REAL)
-                endif
-              endif
-
-              ! source element is poroelastic
-              if (ispec_is_poroelastic(ispec)) then
-!                sourcearray(1,i,j) = real(- sin(anglesource(i_source)) * hlagrange,kind=CUSTOM_REAL)
-!                sourcearray(2,i,j) = real(  cos(anglesource(i_source)) * hlagrange,kind=CUSTOM_REAL)
-!! DK DK May 2018: the sign of the source was inverted compared to the analytical solution for a simple elastic benchmark
-!! DK DK May 2018: with a force source (the example that is in EXAMPLES/check_absolute_amplitude_of_force_source_seismograms),
-!! DK DK May 2018: which means that the sign was not right here. I changed it. Please do NOT revert that change,
-!! DK DK May 2018: otherwise the code will give inverted seismograms compared to analytical solutions for benchmarks,
-!! DK DK May 2018: and more generally compared to reality
-                sourcearray(1,i,j) = real(+ sin(anglesource(i_source)) * hlagrange,kind=CUSTOM_REAL)
-                sourcearray(2,i,j) = real(- cos(anglesource(i_source)) * hlagrange,kind=CUSTOM_REAL)
-              endif
-
-            enddo
-          enddo
-
+          call compute_arrays_source_forcesolution(ispec,hxis,hgammas,sourcearray,anglesource(i_source))
         case (2)
           ! moment-tensor source
-          call compute_arrays_source(ispec,xi_source(i_source),gamma_source(i_source),sourcearray, &
-                                     Mxx(i_source),Mzz(i_source),Mxz(i_source),xix,xiz,gammax,gammaz,xigll,zigll,nspec)
-          ! checks source
-          if (ispec_is_acoustic(ispec)) then
-            call exit_MPI(myrank,'cannot have moment tensor source in acoustic element')
-          endif
-
-          ! checks wave type
-          if (ispec_is_elastic(ispec)) then
-            if (.not. P_SV ) call exit_MPI(myrank,'cannot have moment tensor source in SH (membrane) waves calculation')
-          endif
-
+          call compute_arrays_source_cmt(ispec,hxis,hgammas,hpxis,hpgammas,sourcearray, &
+                                         Mxx(i_source),Mzz(i_source),Mxz(i_source),xix,xiz,gammax,gammaz,nspec)
         end select
 
         ! stores sourcearray for all sources
         sourcearrays(:,:,:,i_source) = sourcearray(:,:,:)
-
       endif
     endif
   enddo