From 47a8f3f250a2cd766c4a41cdb0a11ea2a6b60989 Mon Sep 17 00:00:00 2001 From: Paul Cristini Date: Tue, 10 Nov 2020 14:39:30 +0100 Subject: [PATCH] Update compute_forces_viscoelastic.F90 Correct a bug which broke axisymmetric simulations --- src/specfem2D/compute_forces_viscoelastic.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/specfem2D/compute_forces_viscoelastic.F90 b/src/specfem2D/compute_forces_viscoelastic.F90 index b51ca211b..56a9d4c31 100644 --- a/src/specfem2D/compute_forces_viscoelastic.F90 +++ b/src/specfem2D/compute_forces_viscoelastic.F90 @@ -304,6 +304,7 @@ subroutine compute_forces_viscoelastic(accel_elastic,veloc_elastic,displ_elastic sigma_zz = lambdaplus2mu_unrelaxed_elastic*duz_dzl(i,j) & + lambdal_unrelaxed_elastic * (dux_dxl(i,j) + sigma_zz/xxi) sigma_xz = mul_unrelaxed_elastic*(duz_dxl(i,j) + dux_dzl(i,j)) + sigma_zx = sigma_xz sigma_thetatheta(i,j) = lambdal_unrelaxed_elastic * (duz_dzl(i,j) + dux_dxl(i,j)) & + lambdaplus2mu_unrelaxed_elastic*sigma_thetatheta(i,j)/xxi else @@ -313,6 +314,7 @@ subroutine compute_forces_viscoelastic(accel_elastic,veloc_elastic,displ_elastic sigma_zz = lambdaplus2mu_unrelaxed_elastic*duz_dzl(i,j) & + lambdal_unrelaxed_elastic * (dux_dxl(i,j) + dummy_loc(1,i,j)/coord(1,ibool(i,j,ispec))) sigma_xz = mul_unrelaxed_elastic*(duz_dxl(i,j) + dux_dzl(i,j)) + sigma_zx = sigma_xz sigma_thetatheta(i,j) = lambdal_unrelaxed_elastic * (duz_dzl(i,j) + dux_dxl(i,j)) & + lambdaplus2mu_unrelaxed_elastic * dummy_loc(1,i,j)/coord(1,ibool(i,j,ispec)) r_xiplus1(i,j) = coord(1,ibool(i,j,ispec))/(xiglj(i)+ONE) @@ -324,6 +326,7 @@ subroutine compute_forces_viscoelastic(accel_elastic,veloc_elastic,displ_elastic sigma_zz = lambdaplus2mu_unrelaxed_elastic*duz_dzl(i,j) & + lambdal_unrelaxed_elastic * (dux_dxl(i,j) + dummy_loc(1,i,j)/coord(1,ibool(i,j,ispec))) sigma_xz = mul_unrelaxed_elastic*(duz_dxl(i,j) + dux_dzl(i,j)) + sigma_zx = sigma_xz sigma_thetatheta(i,j) = lambdal_unrelaxed_elastic * (duz_dzl(i,j) + dux_dxl(i,j)) & + lambdaplus2mu_unrelaxed_elastic * dummy_loc(1,i,j)/coord(1,ibool(i,j,ispec)) endif @@ -390,12 +393,14 @@ subroutine compute_forces_viscoelastic(accel_elastic,veloc_elastic,displ_elastic sigma_xx = c11*dux_dxl(i,j) + c13*duz_dzl(i,j) + c12*sigma_xx/xxi sigma_zz = c13*dux_dxl(i,j) + c33*duz_dzl(i,j) + c23*sigma_zz/xxi sigma_xz = c15*dux_dxl(i,j) + c35*duz_dzl(i,j) + c55*(duz_dxl(i,j) + dux_dzl(i,j)) + sigma_zx = sigma_xz sigma_thetatheta(i,j) = c12*dux_dxl(i,j) + c23*duz_dzl(i,j) + c22*sigma_thetatheta(i,j)/xxi else ! not first GLJ point but not on the axis sigma_xx = c11*dux_dxl(i,j) + c13*duz_dzl(i,j) + c12*dummy_loc(1,i,j)/coord(1,ibool(i,j,ispec)) sigma_zz = c13*dux_dxl(i,j) + c33*duz_dzl(i,j) + c23*dummy_loc(1,i,j)/coord(1,ibool(i,j,ispec)) sigma_xz = c15*dux_dxl(i,j) + c35*duz_dzl(i,j) + c55*(duz_dxl(i,j) + dux_dzl(i,j)) + sigma_zx = sigma_xz sigma_thetatheta(i,j) = c12*dux_dxl(i,j) + c23*duz_dzl(i,j) + & c22*dummy_loc(1,i,j)/coord(1,ibool(i,j,ispec)) endif @@ -404,6 +409,7 @@ subroutine compute_forces_viscoelastic(accel_elastic,veloc_elastic,displ_elastic sigma_xx = c11*dux_dxl(i,j) + c13*duz_dzl(i,j) + c12*dummy_loc(1,i,j)/coord(1,ibool(i,j,ispec)) sigma_zz = c13*dux_dxl(i,j) + c33*duz_dzl(i,j) + c23*dummy_loc(1,i,j)/coord(1,ibool(i,j,ispec)) sigma_xz = c15*dux_dxl(i,j) + c35*duz_dzl(i,j) + c55*(duz_dxl(i,j) + dux_dzl(i,j)) + sigma_zx = sigma_xz sigma_thetatheta(i,j) = c12*dux_dxl(i,j) + c23*duz_dzl(i,j) + c22*dummy_loc(1,i,j)/coord(1,ibool(i,j,ispec)) endif else