diff --git a/src/specfem2D/Makefile.in b/src/specfem2D/Makefile.in index e93da9975..40d1b97ba 100644 --- a/src/specfem2D/Makefile.in +++ b/src/specfem2D/Makefile.in @@ -158,6 +158,7 @@ OBJS_SPECFEM2D = \ $O/set_sources.o \ $O/setup_sources_receivers.o \ $O/sort_array_coordinates.o \ + $O/update_displacement_scheme.o \ $O/write_seismograms.o \ $O/specfem2D.o \ $O/write_jpeg_image.o \ @@ -407,6 +408,9 @@ $O/spline_routines.o: ${S}/spline_routines.f90 ${SETUP}/constants.h $O/sort_array_coordinates.o: ${S}/sort_array_coordinates.F90 ${SETUP}/constants.h ${F90} $(FLAGS_CHECK) -c -o $O/sort_array_coordinates.o ${S}/sort_array_coordinates.F90 +$O/update_displacement_scheme.o: ${S}/update_displacement_scheme.F90 ${SETUP}/constants.h + ${F90} $(FLAGS_CHECK) -c -o $O/update_displacement_scheme.o ${S}/update_displacement_scheme.F90 + $O/write_seismograms.o: ${S}/write_seismograms.F90 ${SETUP}/constants.h ${F90} $(FLAGS_CHECK) -c -o $O/write_seismograms.o ${S}/write_seismograms.F90 diff --git a/src/specfem2D/specfem2D.F90 b/src/specfem2D/specfem2D.F90 index c8b2a6674..42c46005c 100644 --- a/src/specfem2D/specfem2D.F90 +++ b/src/specfem2D/specfem2D.F90 @@ -4922,246 +4922,25 @@ program specfem2D do i_stage=1, stage_time_scheme -! update displacement using finite-difference time scheme (Newmark) - if(any_elastic) then - - if(time_stepping_scheme==1)then -#ifdef FORCE_VECTORIZATION -!! DK DK this allows for full vectorization by using a trick to see the 2D array as a 1D array -!! DK DK whose dimension is the product of the two dimensions, the second dimension being equal to 1 - do i = 1,3*nglob_elastic !! DK DK here change 3 to NDIM when/if we suppress the 2nd component of the arrays (the SH component) -! do i = 1,NDIM*nglob_elastic !! DK DK this should be the correct size in principle, but not here because of the SH component - displ_elastic_old(i,1) = displ_elastic(i,1) + deltatsquareover2 * accel_elastic(i,1) - displ_elastic(i,1) = displ_elastic(i,1) & - + deltat*veloc_elastic(i,1) & - + deltatsquareover2*accel_elastic(i,1) - veloc_elastic(i,1) = veloc_elastic(i,1) + deltatover2*accel_elastic(i,1) - enddo -#else - displ_elastic_old = displ_elastic + deltatsquareover2 * accel_elastic - displ_elastic = displ_elastic & - + deltat*veloc_elastic & - + deltatsquareover2*accel_elastic - veloc_elastic = veloc_elastic + deltatover2*accel_elastic -#endif - endif - accel_elastic_adj_coupling = accel_elastic - accel_elastic = ZERO - - if(SIMULATION_TYPE == 3) then ! Adjoint calculation -!! DK DK this should be fully vectorized - b_displ_elastic_old = b_displ_elastic + deltatsquareover2 * b_accel_elastic - b_displ_elastic = b_displ_elastic & - + b_deltat*b_veloc_elastic & - + b_deltatsquareover2*b_accel_elastic - b_veloc_elastic = b_veloc_elastic + b_deltatover2*b_accel_elastic - b_accel_elastic = ZERO - endif - endif - - if(any_poroelastic) then - - if(time_stepping_scheme==1)then - !for the solid - displs_poroelastic_old = displs_poroelastic + deltatover2 * accels_poroelastic - displs_poroelastic = displs_poroelastic & - + deltat*velocs_poroelastic & - + deltatsquareover2*accels_poroelastic - velocs_poroelastic = velocs_poroelastic + deltatover2*accels_poroelastic - accels_poroelastic_adj_coupling = accels_poroelastic - accels_poroelastic = ZERO - !for the fluid - displw_poroelastic = displw_poroelastic & - + deltat*velocw_poroelastic & - + deltatsquareover2*accelw_poroelastic - velocw_poroelastic = velocw_poroelastic + deltatover2*accelw_poroelastic - accelw_poroelastic_adj_coupling = accelw_poroelastic - accelw_poroelastic = ZERO - endif - - if(time_stepping_scheme==2)then - !for the solid - accels_poroelastic_adj_coupling = accels_poroelastic - accels_poroelastic = ZERO - !for the fluid - accelw_poroelastic_adj_coupling = accelw_poroelastic - accelw_poroelastic = ZERO - endif - - if(time_stepping_scheme==3)then - !for the solid - accels_poroelastic_adj_coupling = accels_poroelastic - accels_poroelastic = ZERO - !for the fluid - accelw_poroelastic_adj_coupling = accelw_poroelastic - accelw_poroelastic = ZERO - endif - - if(SIMULATION_TYPE == 3) then ! Adjoint calculation - !for the solid - b_displs_poroelastic = b_displs_poroelastic & - + b_deltat*b_velocs_poroelastic & - + b_deltatsquareover2*b_accels_poroelastic - b_velocs_poroelastic = b_velocs_poroelastic + b_deltatover2*b_accels_poroelastic - b_accels_poroelastic = ZERO - !for the fluid - b_displw_poroelastic = b_displw_poroelastic & - + b_deltat*b_velocw_poroelastic & - + b_deltatsquareover2*b_accelw_poroelastic - b_velocw_poroelastic = b_velocw_poroelastic + b_deltatover2*b_accelw_poroelastic - b_accelw_poroelastic = ZERO - endif - endif - -!-------------------------------------------------------------------------------------------- -! implement viscous attenuation for poroelastic media -! - if(ATTENUATION_PORO_FLUID_PART .and. any_poroelastic) then - ! update memory variables with fourth-order Runge-Kutta time scheme for attenuation - ! loop over spectral elements - do ispec = 1,nspec - - if(poroelastic(ispec) .and. poroelastcoef(2,2,kmato(ispec)) >0.d0) then - - etal_f = poroelastcoef(2,2,kmato(ispec)) - permlxx = permeability(1,kmato(ispec)) - permlxz = permeability(2,kmato(ispec)) - permlzz = permeability(3,kmato(ispec)) - - ! calcul of the inverse of k - - detk = permlxx*permlzz - permlxz*permlxz - - if(detk /= ZERO) then - invpermlxx = permlzz/detk - invpermlxz = -permlxz/detk - invpermlzz = permlxx/detk - else - stop 'Permeability matrix is not invertible' - endif - - ! relaxed viscous coef - bl_unrelaxed_elastic(1) = etal_f*invpermlxx - bl_unrelaxed_elastic(2) = etal_f*invpermlxz - bl_unrelaxed_elastic(3) = etal_f*invpermlzz - - do j=1,NGLLZ - do i=1,NGLLX - - iglob = ibool(i,j,ispec) - viscox_loc(i,j) = velocw_poroelastic(1,iglob)*bl_unrelaxed_elastic(1) + & - velocw_poroelastic(2,iglob) * bl_unrelaxed_elastic(2) - viscoz_loc(i,j) = velocw_poroelastic(1,iglob)*bl_unrelaxed_elastic(2) + & - velocw_poroelastic(2,iglob)*bl_unrelaxed_elastic(3) - - if(time_stepping_scheme == 1) then - ! evolution rx_viscous - Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscox(i,j,ispec) - Snp1 = - (1.d0 - theta_e/theta_s)/theta_s*viscox_loc(i,j) - rx_viscous(i,j,ispec) = alphaval * rx_viscous(i,j,ispec) & - + betaval * Sn + gammaval * Snp1 - - ! evolution rz_viscous - Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscoz(i,j,ispec) - Snp1 = - (1.d0 - theta_e/theta_s)/theta_s*viscoz_loc(i,j) - rz_viscous(i,j,ispec) = alphaval * rz_viscous(i,j,ispec) & - + betaval * Sn + gammaval * Snp1 - endif - - if(time_stepping_scheme == 2) then - Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscox(i,j,ispec) - rx_viscous_LDDRK(i,j,ispec) = alpha_LDDRK(i_stage) * rx_viscous_LDDRK(i,j,ispec) + & - deltat * (Sn + thetainv * rx_viscous(i,j,ispec)) - rx_viscous(i,j,ispec)= rx_viscous(i,j,ispec)+beta_LDDRK(i_stage) * rx_viscous_LDDRK(i,j,ispec) - - Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscoz(i,j,ispec) - rz_viscous_LDDRK(i,j,ispec)= alpha_LDDRK(i_stage) * rz_viscous_LDDRK(i,j,ispec)+& - deltat * (Sn + thetainv * rz_viscous(i,j,ispec)) - rz_viscous(i,j,ispec)= rz_viscous(i,j,ispec)+beta_LDDRK(i_stage) * rz_viscous_LDDRK(i,j,ispec) - endif - - if(time_stepping_scheme == 3) then - - Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscox(i,j,ispec) - rx_viscous_force_RK(i,j,ispec,i_stage) = deltat * (Sn + thetainv * rx_viscous(i,j,ispec)) - - if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then - if(i_stage == 1)weight_rk = 0.5d0 - if(i_stage == 2)weight_rk = 0.5d0 - if(i_stage == 3)weight_rk = 1.0d0 - - if(i_stage==1)then - rx_viscous_initial_rk(i,j,ispec) = rx_viscous(i,j,ispec) - endif - rx_viscous(i,j,ispec) = rx_viscous_initial_rk(i,j,ispec) + & - weight_rk * rx_viscous_force_RK(i,j,ispec,i_stage) - else if(i_stage==4)then - - rx_viscous(i,j,ispec) = rx_viscous_initial_rk(i,j,ispec) + & - 1.0d0 / 6.0d0 * (rx_viscous_force_RK(i,j,ispec,i_stage) + & - 2.0d0 * rx_viscous_force_RK(i,j,ispec,i_stage) + & - 2.0d0 * rx_viscous_force_RK(i,j,ispec,i_stage) + & - rx_viscous_force_RK(i,j,ispec,i_stage)) - endif - - Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscoz(i,j,ispec) - rz_viscous_force_RK(i,j,ispec,i_stage) = deltat * (Sn + thetainv * rz_viscous(i,j,ispec)) - - if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then - if(i_stage == 1)weight_rk = 0.5d0 - if(i_stage == 2)weight_rk = 0.5d0 - if(i_stage == 3)weight_rk = 1.0d0 - if(i_stage==1)then - rz_viscous_initial_rk(i,j,ispec) = rz_viscous(i,j,ispec) - endif - rz_viscous(i,j,ispec) = rz_viscous_initial_rk(i,j,ispec) + & - weight_rk * rz_viscous_force_RK(i,j,ispec,i_stage) - else if(i_stage==4)then - rz_viscous(i,j,ispec) = rz_viscous_initial_rk(i,j,ispec) + & - 1.0d0 / 6.0d0 * (rz_viscous_force_RK(i,j,ispec,i_stage) + & - 2.0d0 * rz_viscous_force_RK(i,j,ispec,i_stage) + & - 2.0d0 * rz_viscous_force_RK(i,j,ispec,i_stage) + & - rz_viscous_force_RK(i,j,ispec,i_stage)) - endif - endif - enddo - enddo - - if(stage_time_scheme == 1) then - ! save visco for Runge-Kutta scheme when used together with Newmark - viscox(:,:,ispec) = viscox_loc(:,:) - viscoz(:,:,ispec) = viscoz_loc(:,:) - endif - - endif ! end of poroelastic element loop + call update_displacement_precondition_newmark(time_stepping_scheme,SIMULATION_TYPE,& + nglob_acoustic,nglob_elastic,nglob_poroelastic,& + any_acoustic,any_elastic,any_poroelastic,deltat,deltatover2,& + deltatsquareover2,potential_acoustic,potential_dot_acoustic,& + potential_dot_dot_acoustic,potential_acoustic_old,& + displ_elastic,veloc_elastic,accel_elastic,displ_elastic_old,& + displs_poroelastic,velocs_poroelastic,accels_poroelastic,& + displs_poroelastic_old,displw_poroelastic,velocw_poroelastic,& + accelw_poroelastic,b_deltat,b_deltatover2,b_deltatsquareover2,& + b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic,& + b_potential_acoustic_old,& + b_displ_elastic,b_veloc_elastic,b_accel_elastic,b_displ_elastic_old,& + accel_elastic_adj_coupling,& + b_displs_poroelastic,b_velocs_poroelastic,b_accels_poroelastic,& + accels_poroelastic_adj_coupling,& + b_displw_poroelastic,b_velocw_poroelastic,b_accelw_poroelastic,& + accelw_poroelastic_adj_coupling) - enddo ! end of spectral element loop - endif ! end of viscous attenuation for porous media - -!----------------------------------------- if(any_acoustic) then - - if(time_stepping_scheme==1)then - ! Newmark time scheme -!! DK DK this should be vectorized - potential_acoustic_old = potential_acoustic + deltatsquareover2*potential_dot_dot_acoustic - potential_acoustic = potential_acoustic + & - deltat*potential_dot_acoustic + deltatsquareover2*potential_dot_dot_acoustic - potential_dot_acoustic = potential_dot_acoustic + deltatover2*potential_dot_dot_acoustic - - endif - potential_dot_dot_acoustic = ZERO - - if(SIMULATION_TYPE == 3) then ! Adjoint calculation -!! DK DK this should be vectorized - b_potential_acoustic_old = b_potential_acoustic + deltatsquareover2*b_potential_dot_dot_acoustic - b_potential_acoustic = b_potential_acoustic + b_deltat*b_potential_dot_acoustic + & - b_deltatsquareover2*b_potential_dot_dot_acoustic - b_potential_dot_acoustic = b_potential_dot_acoustic & - + b_deltatover2*b_potential_dot_dot_acoustic - b_potential_dot_dot_acoustic = ZERO - endif - ! free surface for an acoustic medium if ( nelem_acoustic_surface > 0 ) then call enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, & @@ -6617,6 +6396,133 @@ program specfem2D if(any_poroelastic) then +!-------------------------------------------------------------------------------------------- +! implement viscous attenuation for poroelastic media +! + if(ATTENUATION_PORO_FLUID_PART) then + ! update memory variables with fourth-order Runge-Kutta time scheme for attenuation + ! loop over spectral elements + do ispec = 1,nspec + + if(poroelastic(ispec) .and. poroelastcoef(2,2,kmato(ispec)) >0.d0) then + + etal_f = poroelastcoef(2,2,kmato(ispec)) + permlxx = permeability(1,kmato(ispec)) + permlxz = permeability(2,kmato(ispec)) + permlzz = permeability(3,kmato(ispec)) + + ! calcul of the inverse of k + + detk = permlxx*permlzz - permlxz*permlxz + + if(detk /= ZERO) then + invpermlxx = permlzz/detk + invpermlxz = -permlxz/detk + invpermlzz = permlxx/detk + else + stop 'Permeability matrix is not invertible' + endif + + ! relaxed viscous coef + bl_unrelaxed_elastic(1) = etal_f*invpermlxx + bl_unrelaxed_elastic(2) = etal_f*invpermlxz + bl_unrelaxed_elastic(3) = etal_f*invpermlzz + + do j=1,NGLLZ + do i=1,NGLLX + + iglob = ibool(i,j,ispec) + viscox_loc(i,j) = velocw_poroelastic(1,iglob)*bl_unrelaxed_elastic(1) + & + velocw_poroelastic(2,iglob) * bl_unrelaxed_elastic(2) + viscoz_loc(i,j) = velocw_poroelastic(1,iglob)*bl_unrelaxed_elastic(2) + & + velocw_poroelastic(2,iglob)*bl_unrelaxed_elastic(3) + + if(time_stepping_scheme == 1) then + ! evolution rx_viscous + Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscox(i,j,ispec) + Snp1 = - (1.d0 - theta_e/theta_s)/theta_s*viscox_loc(i,j) + rx_viscous(i,j,ispec) = alphaval * rx_viscous(i,j,ispec) & + + betaval * Sn + gammaval * Snp1 + + ! evolution rz_viscous + Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscoz(i,j,ispec) + Snp1 = - (1.d0 - theta_e/theta_s)/theta_s*viscoz_loc(i,j) + rz_viscous(i,j,ispec) = alphaval * rz_viscous(i,j,ispec) & + + betaval * Sn + gammaval * Snp1 + endif + + if(time_stepping_scheme == 2) then + Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscox(i,j,ispec) + rx_viscous_LDDRK(i,j,ispec) = alpha_LDDRK(i_stage) * rx_viscous_LDDRK(i,j,ispec) + & + deltat * (Sn + thetainv * rx_viscous(i,j,ispec)) + rx_viscous(i,j,ispec)= rx_viscous(i,j,ispec)+beta_LDDRK(i_stage) * rx_viscous_LDDRK(i,j,ispec) + + Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscoz(i,j,ispec) + rz_viscous_LDDRK(i,j,ispec)= alpha_LDDRK(i_stage) * rz_viscous_LDDRK(i,j,ispec)+& + deltat * (Sn + thetainv * rz_viscous(i,j,ispec)) + rz_viscous(i,j,ispec)= rz_viscous(i,j,ispec)+beta_LDDRK(i_stage) * rz_viscous_LDDRK(i,j,ispec) + endif + + if(time_stepping_scheme == 3) then + + Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscox(i,j,ispec) + rx_viscous_force_RK(i,j,ispec,i_stage) = deltat * (Sn + thetainv * rx_viscous(i,j,ispec)) + + if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then + if(i_stage == 1)weight_rk = 0.5d0 + if(i_stage == 2)weight_rk = 0.5d0 + if(i_stage == 3)weight_rk = 1.0d0 + + if(i_stage==1)then + rx_viscous_initial_rk(i,j,ispec) = rx_viscous(i,j,ispec) + endif + rx_viscous(i,j,ispec) = rx_viscous_initial_rk(i,j,ispec) + & + weight_rk * rx_viscous_force_RK(i,j,ispec,i_stage) + else if(i_stage==4)then + + rx_viscous(i,j,ispec) = rx_viscous_initial_rk(i,j,ispec) + & + 1.0d0 / 6.0d0 * (rx_viscous_force_RK(i,j,ispec,i_stage) + & + 2.0d0 * rx_viscous_force_RK(i,j,ispec,i_stage) + & + 2.0d0 * rx_viscous_force_RK(i,j,ispec,i_stage) + & + rx_viscous_force_RK(i,j,ispec,i_stage)) + endif + + Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscoz(i,j,ispec) + rz_viscous_force_RK(i,j,ispec,i_stage) = deltat * (Sn + thetainv * rz_viscous(i,j,ispec)) + + if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then + if(i_stage == 1)weight_rk = 0.5d0 + if(i_stage == 2)weight_rk = 0.5d0 + if(i_stage == 3)weight_rk = 1.0d0 + if(i_stage==1)then + rz_viscous_initial_rk(i,j,ispec) = rz_viscous(i,j,ispec) + endif + rz_viscous(i,j,ispec) = rz_viscous_initial_rk(i,j,ispec) + & + weight_rk * rz_viscous_force_RK(i,j,ispec,i_stage) + else if(i_stage==4)then + rz_viscous(i,j,ispec) = rz_viscous_initial_rk(i,j,ispec) + & + 1.0d0 / 6.0d0 * (rz_viscous_force_RK(i,j,ispec,i_stage) + & + 2.0d0 * rz_viscous_force_RK(i,j,ispec,i_stage) + & + 2.0d0 * rz_viscous_force_RK(i,j,ispec,i_stage) + & + rz_viscous_force_RK(i,j,ispec,i_stage)) + endif + endif + enddo + enddo + + if(stage_time_scheme == 1) then + ! save visco for Runge-Kutta scheme when used together with Newmark + viscox(:,:,ispec) = viscox_loc(:,:) + viscoz(:,:,ispec) = viscoz_loc(:,:) + endif + + endif ! end of poroelastic element loop + + enddo ! end of spectral element loop + endif ! end of viscous attenuation for porous media + +!----------------------------------------- + if(SIMULATION_TYPE == 3) then ! if inviscid fluid, comment the reading and uncomment the zeroing ! read(23,rec=NSTEP-it+1) b_viscodampx diff --git a/src/specfem2D/update_displacement_scheme.F90 b/src/specfem2D/update_displacement_scheme.F90 new file mode 100644 index 000000000..2c250ad32 --- /dev/null +++ b/src/specfem2D/update_displacement_scheme.F90 @@ -0,0 +1,206 @@ +!======================================================================== +! +! S P E C F E M 2 D Version 7 . 0 +! -------------------------------- +! +! Copyright CNRS, INRIA and University of Pau, France, +! and Princeton University / California Institute of Technology, USA. +! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr +! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr +! Roland Martin, roland DOT martin aT univ-pau DOT fr +! Christina Morency, cmorency aT princeton DOT edu +! +! This software is a computer program whose purpose is to solve +! the two-dimensional viscoelastic anisotropic or poroelastic wave equation +! using a spectral-element method (SEM). +! +! This software is governed by the CeCILL license under French law and +! abiding by the rules of distribution of free software. You can use, +! modify and/or redistribute the software under the terms of the CeCILL +! license as circulated by CEA, CNRS and INRIA at the following URL +! "http://www.cecill.info". +! +! As a counterpart to the access to the source code and rights to copy, +! modify and redistribute granted by the license, users are provided only +! with a limited warranty and the software's author, the holder of the +! economic rights, and the successive licensors have only limited +! liability. +! +! In this respect, the user's attention is drawn to the risks associated +! with loading, using, modifying and/or developing or reproducing the +! software by the user in light of its specific status of free software, +! that may mean that it is complicated to manipulate, and that also +! therefore means that it is reserved for developers and experienced +! professionals having in-depth computer knowledge. Users are therefore +! encouraged to load and test the software's suitability as regards their +! requirements in conditions enabling the security of their systems and/or +! data to be ensured and, more generally, to use and operate it in the +! same conditions as regards security. +! +! The full text of the license is available in file "LICENSE". +! +!======================================================================== + + subroutine update_displacement_precondition_newmark(time_stepping_scheme,SIMULATION_TYPE,& + nglob_acoustic,nglob_elastic,nglob_poroelastic,& + any_acoustic,any_elastic,any_poroelastic,deltat,deltatover2,& + deltatsquareover2,potential_acoustic,potential_dot_acoustic,& + potential_dot_dot_acoustic,potential_acoustic_old,& + displ_elastic,veloc_elastic,accel_elastic,displ_elastic_old,& + displs_poroelastic,velocs_poroelastic,accels_poroelastic,& + displs_poroelastic_old,displw_poroelastic,velocw_poroelastic,& + accelw_poroelastic,b_deltat,b_deltatover2,b_deltatsquareover2,& + b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic,& + b_potential_acoustic_old,& + b_displ_elastic,b_veloc_elastic,b_accel_elastic,b_displ_elastic_old,& + accel_elastic_adj_coupling,& + b_displs_poroelastic,b_velocs_poroelastic,b_accels_poroelastic,& + accels_poroelastic_adj_coupling,& + b_displw_poroelastic,b_velocw_poroelastic,b_accelw_poroelastic,& + accelw_poroelastic_adj_coupling) + + implicit none + include 'constants.h' + + integer :: time_stepping_scheme,SIMULATION_TYPE, nglob_acoustic, nglob_elastic,nglob_poroelastic + logical :: any_acoustic, any_elastic, any_poroelastic + +! SIMULATIONTYPE == 1 + double precision :: deltat,deltatover2,deltatsquareover2 + real(kind=CUSTOM_REAL),dimension(nglob_acoustic) :: potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic,& + potential_acoustic_old + real(kind=CUSTOM_REAL),dimension(3,nglob_elastic) :: displ_elastic,veloc_elastic,accel_elastic,displ_elastic_old + real(kind=CUSTOM_REAL),dimension(NDIM,nglob_poroelastic) :: displs_poroelastic,velocs_poroelastic,accels_poroelastic,& + displs_poroelastic_old,& + displw_poroelastic,velocw_poroelastic,accelw_poroelastic + +! SIMULATIONTYPE == 3 + double precision :: b_deltat,b_deltatover2,b_deltatsquareover2 + real(kind=CUSTOM_REAL),dimension(nglob_acoustic) :: b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic,& + b_potential_acoustic_old + real(kind=CUSTOM_REAL),dimension(3,nglob_elastic) :: b_displ_elastic,b_veloc_elastic,b_accel_elastic,b_displ_elastic_old,& + accel_elastic_adj_coupling + real(kind=CUSTOM_REAL),dimension(NDIM,nglob_poroelastic) :: b_displs_poroelastic,b_velocs_poroelastic,b_accels_poroelastic,& + accels_poroelastic_adj_coupling, & + b_displw_poroelastic,b_velocw_poroelastic,b_accelw_poroelastic, & + accelw_poroelastic_adj_coupling + + +! update displacement using finite-difference time scheme (Newmark) + + if(any_acoustic) then + + if(time_stepping_scheme==1)then + ! Newmark time scheme +!! DK DK this should be vectorized + potential_acoustic_old = potential_acoustic + deltatsquareover2*potential_dot_dot_acoustic + potential_acoustic = potential_acoustic + & + deltat*potential_dot_acoustic + deltatsquareover2*potential_dot_dot_acoustic + potential_dot_acoustic = potential_dot_acoustic + deltatover2*potential_dot_dot_acoustic + + endif + potential_dot_dot_acoustic = ZERO + + if(SIMULATION_TYPE == 3) then ! Adjoint calculation +!! DK DK this should be vectorized + b_potential_acoustic_old = b_potential_acoustic + b_deltatsquareover2*b_potential_dot_dot_acoustic + b_potential_acoustic = b_potential_acoustic + b_deltat*b_potential_dot_acoustic + & + b_deltatsquareover2*b_potential_dot_dot_acoustic + b_potential_dot_acoustic = b_potential_dot_acoustic & + + b_deltatover2*b_potential_dot_dot_acoustic + b_potential_dot_dot_acoustic = ZERO + endif + + endif + + if(any_elastic) then + + if(time_stepping_scheme==1)then +#ifdef FORCE_VECTORIZATION +!! DK DK this allows for full vectorization by using a trick to see the 2D array as a 1D array +!! DK DK whose dimension is the product of the two dimensions, the second dimension being equal to 1 + do i = 1,3*nglob_elastic !! DK DK here change 3 to NDIM when/if we suppress the 2nd component of the arrays (the SH component) +! do i = 1,NDIM*nglob_elastic !! DK DK this should be the correct size in principle, but not here because of the SH component + displ_elastic_old(i,1) = displ_elastic(i,1) + deltatsquareover2 * accel_elastic(i,1) + displ_elastic(i,1) = displ_elastic(i,1) & + + deltat*veloc_elastic(i,1) & + + deltatsquareover2*accel_elastic(i,1) + veloc_elastic(i,1) = veloc_elastic(i,1) + deltatover2*accel_elastic(i,1) + enddo +#else + displ_elastic_old = displ_elastic + deltatsquareover2 * accel_elastic + displ_elastic = displ_elastic & + + deltat*veloc_elastic & + + deltatsquareover2*accel_elastic + veloc_elastic = veloc_elastic + deltatover2*accel_elastic +#endif + endif + accel_elastic_adj_coupling = accel_elastic + accel_elastic = ZERO + + if(SIMULATION_TYPE == 3) then ! Adjoint calculation +!! DK DK this should be fully vectorized + b_displ_elastic_old = b_displ_elastic + b_deltatsquareover2 * b_accel_elastic + b_displ_elastic = b_displ_elastic & + + b_deltat*b_veloc_elastic & + + b_deltatsquareover2*b_accel_elastic + b_veloc_elastic = b_veloc_elastic + b_deltatover2*b_accel_elastic + b_accel_elastic = ZERO + endif + endif + + if(any_poroelastic) then + + if(time_stepping_scheme==1)then + !for the solid + displs_poroelastic_old = displs_poroelastic + deltatover2 * accels_poroelastic + displs_poroelastic = displs_poroelastic & + + deltat*velocs_poroelastic & + + deltatsquareover2*accels_poroelastic + velocs_poroelastic = velocs_poroelastic + deltatover2*accels_poroelastic + accels_poroelastic_adj_coupling = accels_poroelastic + accels_poroelastic = ZERO + !for the fluid + displw_poroelastic = displw_poroelastic & + + deltat*velocw_poroelastic & + + deltatsquareover2*accelw_poroelastic + velocw_poroelastic = velocw_poroelastic + deltatover2*accelw_poroelastic + accelw_poroelastic_adj_coupling = accelw_poroelastic + accelw_poroelastic = ZERO + endif + + if(time_stepping_scheme==2)then + !for the solid + accels_poroelastic_adj_coupling = accels_poroelastic + accels_poroelastic = ZERO + !for the fluid + accelw_poroelastic_adj_coupling = accelw_poroelastic + accelw_poroelastic = ZERO + endif + + if(time_stepping_scheme==3)then + !for the solid + accels_poroelastic_adj_coupling = accels_poroelastic + accels_poroelastic = ZERO + !for the fluid + accelw_poroelastic_adj_coupling = accelw_poroelastic + accelw_poroelastic = ZERO + endif + + if(SIMULATION_TYPE == 3) then ! Adjoint calculation + !for the solid + b_displs_poroelastic = b_displs_poroelastic & + + b_deltat*b_velocs_poroelastic & + + b_deltatsquareover2*b_accels_poroelastic + b_velocs_poroelastic = b_velocs_poroelastic + b_deltatover2*b_accels_poroelastic + b_accels_poroelastic = ZERO + !for the fluid + b_displw_poroelastic = b_displw_poroelastic & + + b_deltat*b_velocw_poroelastic & + + b_deltatsquareover2*b_accelw_poroelastic + b_velocw_poroelastic = b_velocw_poroelastic + b_deltatover2*b_accelw_poroelastic + b_accelw_poroelastic = ZERO + endif + endif + + end subroutine update_displacement_precondition_newmark