[cig-commits] r21768 - seismo/3D/SPECFEM3D/branches/coupling_vadim/src/specfem3D

[vadim at geodynamics.org]

Author: vadim
Date: 2013-04-08 10:24:52 -0700 (Mon, 08 Apr 2013)
New Revision: 21768

Modified:
   seismo/3D/SPECFEM3D/branches/coupling_vadim/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
   seismo/3D/SPECFEM3D/branches/coupling_vadim/src/specfem3D/compute_stacey_viscoelastic.f90
Log:
add dsm traction reading in viscoelastic case


Modified: seismo/3D/SPECFEM3D/branches/coupling_vadim/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D/branches/coupling_vadim/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90	2013-04-08 17:19:04 UTC (rev 21767)
+++ seismo/3D/SPECFEM3D/branches/coupling_vadim/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90	2013-04-08 17:24:52 UTC (rev 21768)
@@ -240,7 +240,7 @@
                         coupling_el_po_jacobian2Dw, &
                         ispec_is_inner,phase_is_inner)
     endif
-
+    if (not.USE_VADIM) then
 ! adds source term (single-force/moment-tensor solution)
     call compute_add_sources_viscoelastic( NSPEC_AB,NGLOB_AB,accel, &
                         ibool,ispec_is_inner,phase_is_inner, &
@@ -251,7 +251,7 @@
                         nadj_rec_local,adj_sourcearrays,b_accel, &
                         NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
                         GPU_MODE, Mesh_pointer )
-
+    end if
     ! assemble all the contributions between slices using MPI
     if( phase_is_inner .eqv. .false. ) then
        ! sends accel values to corresponding MPI interface neighbors

Modified: seismo/3D/SPECFEM3D/branches/coupling_vadim/src/specfem3D/compute_stacey_viscoelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/coupling_vadim/src/specfem3D/compute_stacey_viscoelastic.f90	2013-04-08 17:19:04 UTC (rev 21767)
+++ seismo/3D/SPECFEM3D/branches/coupling_vadim/src/specfem3D/compute_stacey_viscoelastic.f90	2013-04-08 17:24:52 UTC (rev 21768)
@@ -84,6 +84,21 @@
   integer :: ispec,iglob,i,j,k,iface,igll
   !integer:: reclen1,reclen2
 
+! VM VM for DSM/SEM Hybrid method
+  integer Ntime_step_dsm
+  real(kind=CUSTOM_REAL) :: Veloc_dsm_boundary(3,Ntime_step_dsm,NGLLSQUARE,num_abs_boundary_faces)
+  real(kind=CUSTOM_REAL) :: Tract_dsm_boundary(3,Ntime_step_dsm,NGLLSQUARE,num_abs_boundary_faces)
+
+  if (USE_VADIM) then
+     if ( phase_is_inner .eqv. .false. ) then
+        if (mod(it_dsm,Ntime_step_dsm+1)==0.or.it==1) then
+           call read_dsm_file(Veloc_dsm_boundary,Tract_dsm_boundary,&
+                Ntime_step_dsm,num_abs_boundary_faces,&
+                IIN_veloc_dsm,IIN_tract_dsm,it_dsm)
+        end if
+     end if
+  end if
+
   ! checks if anything to do
   if( num_abs_boundary_faces == 0 ) return
 
@@ -126,7 +141,11 @@
                  vx=veloc(1,iglob)
                  vy=veloc(2,iglob)
                  vz=veloc(3,iglob)
-
+                 if (USE_VADIM) then
+                     vx = vx - Veloc_dsm_boundary(1,it_dsm,igll,iface)
+                     vy = vy - Veloc_dsm_boundary(2,it_dsm,igll,iface)
+                     vz = vz - Veloc_dsm_boundary(3,it_dsm,igll,iface)
+                 end if
                  ! gets associated normal
                  nx = abs_boundary_normal(1,igll,iface)
                  ny = abs_boundary_normal(2,igll,iface)
@@ -140,6 +159,12 @@
                  ty = rho_vp(i,j,k,ispec)*vn*ny + rho_vs(i,j,k,ispec)*(vy-vn*ny)
                  tz = rho_vp(i,j,k,ispec)*vn*nz + rho_vs(i,j,k,ispec)*(vz-vn*nz)
 
+                 if (USE_VADIM) then
+                     tx = -Tract_dsm_boundary(1,it_dsm,igll,iface) + tx 
+                     ty = -Tract_dsm_boundary(2,it_dsm,igll,iface) + ty
+                     tz = -Tract_dsm_boundary(3,it_dsm,igll,iface) + tz 
+                 end if
+
                  ! gets associated, weighted jacobian
                  jacobianw = abs_boundary_jacobian2Dw(igll,iface)
 
@@ -182,3 +207,35 @@
 
   end subroutine compute_stacey_viscoelastic
 
+
+  subroutine read_dsm_file(Veloc_dsm_boundary,Tract_dsm_boundary,&
+         Ntime_step_dsm,num_abs_boundary_faces,&
+         IIN_veloc_dsm,IIN_tract_dsm,it_dsm)
+
+
+   implicit none
+   include "constants.h"
+   integer Ntime_step_dsm,igll,it_dsm
+   integer iface,num_abs_boundary_faces,i,j,IIN_veloc_dsm,IIN_tract_dsm
+   real(kind=CUSTOM_REAL) :: Veloc_dsm_boundary(3,Ntime_step_dsm,NGLLSQUARE,num_abs_boundary_faces)
+   real(kind=CUSTOM_REAL) :: Tract_dsm_boundary(3,Ntime_step_dsm,NGLLSQUARE,num_abs_boundary_faces)
+   real(kind=CUSTOM_REAL) :: dsm_boundary_tmp(3,100,5,5)  !!! warning code en dur
+
+   it_dsm = 1
+   write(*,*) 'read dsm files',it_dsm
+   do iface=1,num_abs_boundary_faces
+      
+      igll = 0
+      do j=1,5
+        do i=1,5
+           igll = igll + 1
+           read(IIN_veloc_dsm) dsm_boundary_tmp(:,:,i,j)
+           Veloc_dsm_boundary(:,:,igll,iface) = dsm_boundary_tmp(:,:,i,j) 
+           read(IIN_tract_dsm) dsm_boundary_tmp(:,:,i,j)
+           Tract_dsm_boundary(:,:,igll,iface) = dsm_boundary_tmp(:,:,i,j) 
+      end do
+    end do
+   end do
+
+  
+  end subroutine read_dsm_file