[cig-commits] r19871 - in seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src: cuda meshfem3D shared specfem3D

danielpeter at geodynamics.org danielpeter at geodynamics.org
Mon Mar 26 08:10:51 PDT 2012


Author: danielpeter
Date: 2012-03-26 08:10:51 -0700 (Mon, 26 Mar 2012)
New Revision: 19871

Modified:
   seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
   seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/save_and_compare_cpu_vs_gpu.c
   seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/define_subregions.f90
   seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/constants.h.in
   seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_acoustic.f90
   seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_elastic.f90
   seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_interpolated_dva.f90
   seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_kernels.f90
   seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90
Log:
adds flag USE_RICKER_IPATI to use ricker source time function as default rather than a gaussian

Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu	2012-03-26 14:43:02 UTC (rev 19870)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu	2012-03-26 15:10:51 UTC (rev 19871)
@@ -102,6 +102,9 @@
       fprintf(stderr,"Error after %s: %s\n", kernel_name, cudaGetErrorString(err));
       pause_for_debugger(0);
       //free(kernel_name);
+#ifdef WITH_MPI
+      MPI_Abort(MPI_COMM_WORLD,1);
+#endif      
       exit(EXIT_FAILURE);
     }
 }

Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/save_and_compare_cpu_vs_gpu.c
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/save_and_compare_cpu_vs_gpu.c	2012-03-26 14:43:02 UTC (rev 19870)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/save_and_compare_cpu_vs_gpu.c	2012-03-26 15:10:51 UTC (rev 19871)
@@ -37,6 +37,7 @@
 
 #define MAX(a, b) (((a) > (b)) ? (a) : (b))
 
+
 void save_to_max_surface_file_(float* maxval) {
   int rank;
   char filename[BUFSIZ];
@@ -52,6 +53,7 @@
   fclose(fp);
 }
 
+
 void save_fvector_(float* vector, int* size, int* id, int* cpu_or_gpu) {
   FILE* fp;
   char filename[BUFSIZ];
@@ -265,7 +267,7 @@
   for(i=0;i<*size;i++) {
     if((fabs(vector[i] - compare_vector[i])/vector[i] > 0.0001)) {
       if(error_count < 30) {
-  printf("ERROR[%d]: %g != %g\n",i,compare_vector[i], vector[i]);
+        printf("ERROR[%d]: %f != %f\n",i,compare_vector[i], vector[i]);
       }
       error_count++;
       /* if(compare_vector[i] > 1e-30) error_count++; */
@@ -308,7 +310,7 @@
   int error_count=0;
   for(i=0;i<*size;i++) {
     if((abs(vector[i] - compare_vector[i])/vector[i] > 0.01) && error_count < 30) {
-      printf("ERROR[%d]: %g != %g\n",i,compare_vector[i], vector[i]);
+      printf("ERROR[%d]: %d != %d\n",i,compare_vector[i], vector[i]);
       error_count++;
     }
   }

Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/define_subregions.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/define_subregions.f90	2012-03-26 14:43:02 UTC (rev 19870)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/define_subregions.f90	2012-03-26 15:10:51 UTC (rev 19871)
@@ -72,7 +72,7 @@
      iy1=2*(subregions(isubregion,3) - iproc_eta*NEX_PER_PROC_ETA - 1)
      if(iy1 < 0) iy1 = 0
      iy2=2*(subregions(isubregion,4) - iproc_eta*NEX_PER_PROC_ETA - 1)
-     if(iy2 > 2*(NEX_PER_PROC_XI - 1)) iy2 = 2*(NEX_PER_PROC_ETA - 1)
+     if(iy2 > 2*(NEX_PER_PROC_ETA - 1)) iy2 = 2*(NEX_PER_PROC_ETA - 1)
      diy=2
 
      ir1=2*(subregions(isubregion,5) - 1)

Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/constants.h.in	2012-03-26 14:43:02 UTC (rev 19870)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/constants.h.in	2012-03-26 15:10:51 UTC (rev 19871)
@@ -199,6 +199,9 @@
   double precision, parameter :: FACTOR_FORCE_SOURCE = 1.d15
   integer, parameter :: COMPONENT_FORCE_SOURCE = 3  ! takes direction in comp E/N/Z = 1/2/3
 
+! set to use a Ricker source time function instead of a gaussian
+  logical, parameter :: USE_RICKER_IPATI = .false.
+
 ! use this t0 as earliest starting time rather than the automatically calculated one
 ! (must be positive and bigger than the automatically one to be effective;
 !  simulation will start at t = - t0)

Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_acoustic.f90	2012-03-26 14:43:02 UTC (rev 19870)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_acoustic.f90	2012-03-26 15:10:51 UTC (rev 19871)
@@ -83,7 +83,8 @@
   real(kind=CUSTOM_REAL),dimension(NGLOB_ADJOINT):: b_potential_dot_dot_acoustic
   logical :: ibool_read_adj_arrays
   integer :: it_sub_adj,itime,NTSTEP_BETWEEN_READ_ADJSRC
-  real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
+  real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: &
+    adj_sourcearrays
 
 ! local parameters
   double precision :: f0
@@ -124,8 +125,13 @@
             stf_pre_compute(isource) = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr( &
                   dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
           else
-            stf_pre_compute(isource) = comp_source_time_function_gauss( &
-                  dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+            if( USE_RICKER_IPATI ) then
+              stf_pre_compute(isource) = comp_source_time_function_rickr( &
+                          dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
+            else
+              stf_pre_compute(isource) = comp_source_time_function_gauss( &
+                          dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+            endif
           endif
         enddo
         stf_used_total = stf_used_total + sum(stf_pre_compute(:))
@@ -172,7 +178,8 @@
 
                 ! we use nu_source(:,3) here because we want a source normal to the surface.
                 ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
-                stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+                stf_used = FACTOR_FORCE_SOURCE * &
+                           comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
 
                 ! 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
@@ -187,8 +194,14 @@
 
               else
 
-                ! gaussian source time
-                stf = comp_source_time_function_gauss(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+                if( USE_RICKER_IPATI ) then
+                  stf = comp_source_time_function_rickr( &
+                                dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
+                else
+                  ! gaussian source time
+                  stf = comp_source_time_function_gauss( &
+                                dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+                endif
 
                 ! quasi-heaviside
                 !stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
@@ -306,8 +319,9 @@
            it_start = NSTEP - it_sub_adj*NTSTEP_BETWEEN_READ_ADJSRC + 1
            it_end   = it_start + NTSTEP_BETWEEN_READ_ADJSRC - 1
            write(procname,"(i4)") myrank
-           open(unit=IIN_SU1, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj', &
-                              access='direct',recl=240+4*(NSTEP))
+           open(unit=IIN_SU1, &
+                file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj', &
+                access='direct',recl=240+4*(NSTEP))
            do irec_local = 1,nrec_local
              irec = number_receiver_global(irec_local)
              read(IIN_SU1,rec=irec_local) r4head, adj_temp
@@ -403,8 +417,13 @@
             stf_pre_compute(isource) = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr( &
                   dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
           else
-            stf_pre_compute(isource) = comp_source_time_function_gauss( &
-                  dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+            if( USE_RICKER_IPATI ) then
+              stf_pre_compute(isource) = comp_source_time_function_rickr( &
+                            dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
+            else
+              stf_pre_compute(isource) = comp_source_time_function_gauss( &
+                            dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+            endif
           endif
         enddo
         stf_used_total = stf_used_total + sum(stf_pre_compute(:))
@@ -466,8 +485,14 @@
 
               else
 
-                ! gaussian source time
-                stf = comp_source_time_function_gauss(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+                if( USE_RICKER_IPATI ) then
+                  stf = comp_source_time_function_rickr( &
+                                  dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
+                else
+                  ! gaussian source time
+                  stf = comp_source_time_function_gauss( &
+                                  dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+                endif
 
                 ! distinguishes between single and double precision for reals
                 if(CUSTOM_REAL == SIZE_REAL) then

Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_elastic.f90	2012-03-26 14:43:02 UTC (rev 19870)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_elastic.f90	2012-03-26 15:10:51 UTC (rev 19871)
@@ -91,7 +91,8 @@
   real(kind=CUSTOM_REAL),dimension(NDIM,NGLOB_ADJOINT):: b_accel
   logical :: ibool_read_adj_arrays
   integer :: it_sub_adj,itime,NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY
-  real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
+  real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: &
+    adj_sourcearrays
 
 ! local parameters
   double precision :: f0
@@ -126,8 +127,13 @@
 
     if(GPU_MODE) then
       do isource = 1,NSOURCES
-        stf_pre_compute(isource) = comp_source_time_function( &
+        if( USE_RICKER_IPATI ) then
+          stf_pre_compute(isource) = comp_source_time_function_rickr( &
+                                        dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))          
+        else
+          stf_pre_compute(isource) = comp_source_time_function( &
                                         dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+        endif
       enddo
       ! only implements SIMTYPE=1 and NOISE_TOM=0
       ! write(*,*) "fortran dt = ", dt
@@ -151,49 +157,54 @@
 
                   if(USE_FORCE_POINT_SOURCE) then
 
-                     ! note: for use_force_point_source xi/eta/gamma are in the range [1,NGLL*]
-                     iglob = ibool(nint(xi_source(isource)), &
+                    ! note: for use_force_point_source xi/eta/gamma are in the range [1,NGLL*]
+                    iglob = ibool(nint(xi_source(isource)), &
                           nint(eta_source(isource)), &
                           nint(gamma_source(isource)), &
                           ispec_selected_source(isource))
 
-                     f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
+                    f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
 
-                     !if (it == 1 .and. myrank == 0) then
-                     !  write(IMAIN,*) 'using a source of dominant frequency ',f0
-                     !  write(IMAIN,*) 'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
-                     !  write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
-                     !endif
+                    !if (it == 1 .and. myrank == 0) then
+                    !  write(IMAIN,*) 'using a source of dominant frequency ',f0
+                    !  write(IMAIN,*) 'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+                    !  write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+                    !endif
 
-                     ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
-                     stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+                    ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
+                    stf_used = FACTOR_FORCE_SOURCE * &
+                               comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
 
-                     ! we use a force in a single direction along one of the components:
-                     !  x/y/z or E/N/Z-direction would correspond to 1/2/3 = COMPONENT_FORCE_SOURCE
-                     ! e.g. nu_source(:,3) here would be a source normal to the surface (z-direction).
-                     accel(:,iglob) = accel(:,iglob)  &
+                    ! we use a force in a single direction along one of the components:
+                    !  x/y/z or E/N/Z-direction would correspond to 1/2/3 = COMPONENT_FORCE_SOURCE
+                    ! e.g. nu_source(:,3) here would be a source normal to the surface (z-direction).
+                    accel(:,iglob) = accel(:,iglob)  &
                           + sngl( nu_source(COMPONENT_FORCE_SOURCE,:,isource) ) * stf_used
 
                   else
 
-                     stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+                    if( USE_RICKER_IPATI) then
+                      stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
+                    else
+                      stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+                    endif
+                     
+                    !     distinguish between single and double precision for reals
+                    if(CUSTOM_REAL == SIZE_REAL) then
+                      stf_used = sngl(stf)
+                    else
+                      stf_used = stf
+                    endif
 
-                     !     distinguish between single and double precision for reals
-                     if(CUSTOM_REAL == SIZE_REAL) then
-                        stf_used = sngl(stf)
-                     else
-                        stf_used = stf
-                     endif
-
-                     !     add source array
-                     do k=1,NGLLZ
-                        do j=1,NGLLY
-                           do i=1,NGLLX
-                              iglob = ibool(i,j,k,ispec)
-                              accel(:,iglob) = accel(:,iglob) + sourcearrays(isource,:,i,j,k)*stf_used
-                           enddo
+                    !     add source array
+                    do k=1,NGLLZ
+                      do j=1,NGLLY
+                        do i=1,NGLLX
+                          iglob = ibool(i,j,k,ispec)
+                          accel(:,iglob) = accel(:,iglob) + sourcearrays(isource,:,i,j,k)*stf_used
                         enddo
-                     enddo
+                      enddo
+                    enddo
 
                   endif ! USE_FORCE_POINT_SOURCE
 
@@ -379,8 +390,13 @@
 
     if(GPU_MODE) then
       do isource = 1,NSOURCES
-        stf_pre_compute(isource) = comp_source_time_function( &
-                                          dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+        if( USE_RICKER_IPATI ) then
+          stf_pre_compute(isource) = comp_source_time_function_rickr( &
+                                      dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
+        else
+          stf_pre_compute(isource) = comp_source_time_function( &
+                                      dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+        endif
       enddo
 
       call compute_add_sources_el_s3_cuda(Mesh_pointer,stf_pre_compute, &
@@ -417,7 +433,8 @@
                      !endif
 
                      ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
-                     stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
+                     stf_used = FACTOR_FORCE_SOURCE * &
+                                comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
 
                      ! e.g. we use nu_source(:,3) here if we want a source normal to the surface.
                      ! note: time step is now at NSTEP-it
@@ -428,8 +445,14 @@
 
                      ! see note above: time step corresponds now to NSTEP-it
                      ! (also compare to it-1 for forward simulation)
-                     stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
-
+                     if( USE_RICKER_IPATI ) then
+                       stf = comp_source_time_function_rickr( &
+                                      dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
+                     else
+                       stf = comp_source_time_function( &
+                                      dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+                     endif
+                     
                      ! distinguish between single and double precision for reals
                      if(CUSTOM_REAL == SIZE_REAL) then
                         stf_used = sngl(stf)

Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_interpolated_dva.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_interpolated_dva.f90	2012-03-26 14:43:02 UTC (rev 19870)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_interpolated_dva.f90	2012-03-26 15:10:51 UTC (rev 19871)
@@ -203,8 +203,8 @@
           dzd = dzd + hlagrange*displ_element(3,i,j,k)
           ! velocity
           vxd = vxd + hlagrange*veloc_element(1,i,j,k)
-          vyd = vxd + hlagrange*veloc_element(2,i,j,k)
-          vzd = vxd + hlagrange*veloc_element(3,i,j,k)
+          vyd = vyd + hlagrange*veloc_element(2,i,j,k)
+          vzd = vzd + hlagrange*veloc_element(3,i,j,k)
 
           ! x component -> acoustic potential
           axd = axd + hlagrange*potential_acoustic(iglob)

Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_kernels.f90	2012-03-26 14:43:02 UTC (rev 19870)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_kernels.f90	2012-03-26 15:10:51 UTC (rev 19871)
@@ -197,8 +197,9 @@
                       xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
                       ibool,rhostore,GRAVITY)
       ! adjoint fields: acceleration vector
+      ! new expression (\partial_t^2\bfs^\dagger=-\frac{1}{\rho}\bfnabla\phi^\dagger)
       call compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
-                      potential_dot_dot_acoustic, accel_elm,&
+                      potential_acoustic, accel_elm,&
                       hprime_xx,hprime_yy,hprime_zz, &
                       xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
                       ibool,rhostore,GRAVITY)
@@ -207,17 +208,18 @@
         do j = 1, NGLLY
           do i = 1, NGLLX
             iglob = ibool(i,j,k,ispec)
-
+            
+            ! new expression
             ! density kernel
             rhol = rhostore(i,j,k,ispec)
             rho_ac_kl(i,j,k,ispec) =  rho_ac_kl(i,j,k,ispec) &
-                      - deltat * rhol * dot_product(accel_elm(:,i,j,k), b_displ_elm(:,i,j,k))
+                      + deltat * rhol * dot_product(accel_elm(:,i,j,k), b_displ_elm(:,i,j,k))
 
             ! bulk modulus kernel
-            kappal = kappastore(i,j,k,ispec)
+            kappal = 1._CUSTOM_REAL / kappastore(i,j,k,ispec)
             kappa_ac_kl(i,j,k,ispec) = kappa_ac_kl(i,j,k,ispec) &
-                                  - deltat / kappal  &
-                                  * potential_dot_dot_acoustic(iglob) &
+                                  + deltat * kappal  &
+                                  * potential_acoustic(iglob) &
                                   * b_potential_dot_dot_acoustic(iglob)
 
           enddo

Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90	2012-03-26 14:43:02 UTC (rev 19870)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90	2012-03-26 15:10:51 UTC (rev 19871)
@@ -129,7 +129,7 @@
   ! note: an earlier start time also reduces numerical noise due to a
   !          non-zero offset at the beginning of the source time function
   t0 = - 2.0d0 * minval(tshift_cmt(:) - hdur(:))   ! - 1.5d0 * minval(tshift_cmt-hdur)
-
+  
   ! uses an earlier start time if source is acoustic with a gaussian source time function
   t0_acoustic = 0.0d0
   do isource = 1,NSOURCES
@@ -148,7 +148,7 @@
   call max_all_all_dp(t0_acoustic,t0)
 
   ! point force sources will start depending on the frequency given by hdur
-  if( USE_FORCE_POINT_SOURCE ) then
+  if( USE_FORCE_POINT_SOURCE .or. USE_RICKER_IPATI ) then
     ! note: point force sources will give the dominant frequency in hdur,
     !          thus the main period is 1/hdur.
     !          also, these sources use a Ricker source time function instead of a gaussian.



More information about the CIG-COMMITS mailing list