[cig-commits] [commit] devel: fixed the attenuation bug detected by Paul Cristini: attenuation factors were computed in acoustic (fluid) layers even though attenuation is not implemented there, resulting in meaningless attenuation parameters that in turn could make the code become unstable in some cases (38da483)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Sun Jul 13 18:58:50 PDT 2014
Repository : https://github.com/geodynamics/specfem2d
On branch : devel
Link : https://github.com/geodynamics/specfem2d/compare/4ee2f13012a012f795f6e9ba40582b999af77581...38da48351ad4c6f04247e1cba0937d9af1861b43
>---------------------------------------------------------------
commit 38da48351ad4c6f04247e1cba0937d9af1861b43
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Mon Jul 14 03:51:04 2014 +0200
fixed the attenuation bug detected by Paul Cristini: attenuation factors were computed in acoustic (fluid) layers even though attenuation is not implemented there, resulting in meaningless attenuation parameters that in turn could make the code become unstable in some cases
>---------------------------------------------------------------
38da48351ad4c6f04247e1cba0937d9af1861b43
setup/constants.h.in | 2 +-
src/meshfem2D/meshfem2D.F90 | 3 +++
src/meshfem2D/read_interfaces_file.f90 | 1 -
src/specfem2D/attenuation_model.f90 | 17 +++++++++++++++-
src/specfem2D/compute_forces_viscoelastic.F90 | 17 +++++++++++++---
src/specfem2D/read_external_model.f90 | 28 ++++++++++++++++++++++----
src/specfem2D/specfem2D.F90 | 29 +++++++++++++++++++++++----
7 files changed, 83 insertions(+), 14 deletions(-)
diff --git a/setup/constants.h.in b/setup/constants.h.in
index 3f0cfdf..fb4bff8 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -40,7 +40,7 @@
! For more details, see the comment about this at the beginning of the "SolvOpt" routine
! in file src/specfem2D/attenuation_model.f90.
! We thus strongly discourage using the old routine.
- logical, parameter :: USE_NEW_SOLVOPT_ROUTINE = .false. !! .true.
+ logical, parameter :: USE_NEW_SOLVOPT_ROUTINE = .true.
! use SolvOpt nonlinear optimization with constraints or classical linear least squares
logical, parameter :: USE_SOLVOPT_INSTEAD_OF_LINEAR = .true.
diff --git a/src/meshfem2D/meshfem2D.F90 b/src/meshfem2D/meshfem2D.F90
index a9d0435..75c61e5 100644
--- a/src/meshfem2D/meshfem2D.F90
+++ b/src/meshfem2D/meshfem2D.F90
@@ -1090,6 +1090,9 @@ program meshfem2D
endif
print *
+ if(associated(nz_layer)) deallocate(nz_layer)
+ if(associated(elmnts)) deallocate(elmnts)
+
end program meshfem2D
diff --git a/src/meshfem2D/read_interfaces_file.f90 b/src/meshfem2D/read_interfaces_file.f90
index 36f335b..170ee32 100644
--- a/src/meshfem2D/read_interfaces_file.f90
+++ b/src/meshfem2D/read_interfaces_file.f90
@@ -174,7 +174,6 @@ contains
endif
-
end subroutine read_interfaces_file
end module interfaces_file
diff --git a/src/specfem2D/attenuation_model.f90 b/src/specfem2D/attenuation_model.f90
index 38a3be0..94dad6e 100644
--- a/src/specfem2D/attenuation_model.f90
+++ b/src/specfem2D/attenuation_model.f90
@@ -835,6 +835,7 @@ SUBROUTINE solvopt(n,x,f,fun,flg,grad,options,flfc,func,flgc,gradc,Qref,Kopt,the
if (n<2) then
print *, 'SolvOpt error:'
print *, 'Improper space dimension.'
+ stop 'error in allocate statement in SolvOpt'
options(9)=-one
goto 999
endif
@@ -844,71 +845,85 @@ SUBROUTINE solvopt(n,x,f,fun,flg,grad,options,flfc,func,flgc,gradc,Qref,Kopt,the
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (g(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (g0(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (g1(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (gt(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (gc(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (z(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (x1(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (xopt(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (xrec(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (grec(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (xx(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (deltax(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
allocate (idx(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
+ stop 'error in allocate statement in SolvOpt'
endif
! store flags:
@@ -1943,7 +1958,7 @@ SUBROUTINE solvopt(n,x,f,fun,flg,grad,options,flfc,func,flgc,gradc,Qref,Kopt,the
999 continue
! deallocate working arrays:
- deallocate (idx,deltax,xx,grec,xrec,xopt,x1,z,gc,gt,g1,g0,g,B)
+ deallocate (idx,deltax,xx,grec,xrec,xopt,x1,z,gc,gt,g1,g0,g,B)
END SUBROUTINE solvopt
diff --git a/src/specfem2D/compute_forces_viscoelastic.F90 b/src/specfem2D/compute_forces_viscoelastic.F90
index 584a7b1..9d96535 100644
--- a/src/specfem2D/compute_forces_viscoelastic.F90
+++ b/src/specfem2D/compute_forces_viscoelastic.F90
@@ -67,7 +67,7 @@ subroutine compute_forces_viscoelastic(p_sv,nglob,nspec,myrank,nelemabs,numat, &
rmemory_displ_elastic,rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz, &
rmemory_dux_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dx_prime,rmemory_duz_dz_prime, &
rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
- PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,backward_simulation,STACEY_BOUNDARY_CONDITIONS)
+ PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,backward_simulation,STACEY_BOUNDARY_CONDITIONS,acoustic)
! compute forces for the elastic elements
@@ -83,6 +83,7 @@ subroutine compute_forces_viscoelastic(p_sv,nglob,nspec,myrank,nelemabs,numat, &
integer :: nrec,SIMULATION_TYPE
integer, dimension(nrec) :: ispec_selected_rec,which_proc_receiver
+ logical, dimension(nspec) :: acoustic
integer :: nspec_left,nspec_right,nspec_bottom,nspec_top
integer, dimension(nelemabs) :: ib_left
integer, dimension(nelemabs) :: ib_right
@@ -263,8 +264,17 @@ subroutine compute_forces_viscoelastic(p_sv,nglob,nspec,myrank,nelemabs,numat, &
! loop over spectral elements
do ispec = 1,nspec
+
+! attenuation is not implemented in acoustic (i.e. fluid) media for now, only in viscoelastic (i.e. solid) media
+ if(acoustic(ispec)) cycle
+
if((.not. PML_BOUNDARY_CONDITIONS) .or. (PML_BOUNDARY_CONDITIONS .and. (.not. is_PML(ispec))))then
- do j=1,NGLLZ; do i=1,NGLLX
+ do j=1,NGLLZ
+ do i=1,NGLLX
+
+ ! convention to indicate that Q = 9999 in that element i.e. that there is no viscoelasticity in that element
+ if(inv_tau_sigma_nu1(i,j,ispec,1) < 0.) cycle
+
theta_n_u = dux_dxl_n(i,j,ispec) + duz_dzl_n(i,j,ispec)
theta_nsub1_u = dux_dxl_nsub1(i,j,ispec) + duz_dzl_nsub1(i,j,ispec)
@@ -363,7 +373,8 @@ subroutine compute_forces_viscoelastic(p_sv,nglob,nspec,myrank,nelemabs,numat, &
endif
endif
enddo
- enddo; enddo
+ enddo
+ enddo
endif
enddo
endif
diff --git a/src/specfem2D/read_external_model.f90 b/src/specfem2D/read_external_model.f90
index 699718b..f89d0dc 100644
--- a/src/specfem2D/read_external_model.f90
+++ b/src/specfem2D/read_external_model.f90
@@ -106,8 +106,8 @@
! vsext(i,j,ispec)=0.0
! QKappa, Qmu : dummy values. If attenuation needed than the "read" line and model_velocity.dat_input
! need to be modified to provide QKappa & Qmu values
- QKappa_attenuationext(i,j,ispec) = 10.d0
- Qmu_attenuationext(i,j,ispec) = 10.d0
+ QKappa_attenuationext(i,j,ispec) = 9999.d0
+ Qmu_attenuationext(i,j,ispec) = 9999.d0
enddo
enddo
enddo
@@ -160,6 +160,15 @@
elastic(:) = .false.
poroelastic(:) = .false.
+! initialize to dummy values
+! convention to indicate that Q = 9999 in that element i.e. that there is no viscoelasticity in that element
+ inv_tau_sigma_nu1(:,:,:,:) = -1._CUSTOM_REAL
+ phi_nu1(:,:,:,:) = -1._CUSTOM_REAL
+ inv_tau_sigma_nu2(:,:,:,:) = -1._CUSTOM_REAL
+ phi_nu2(:,:,:,:) = -1._CUSTOM_REAL
+ Mu_nu1(:,:,:) = -1._CUSTOM_REAL
+ Mu_nu2(:,:,:) = -1._CUSTOM_REAL
+
do ispec = 1,nspec
previous_vsext = -1.d0
@@ -178,8 +187,8 @@
poroelastic(ispec) = .false.
elastic(ispec) = .true.
any_elastic = .true.
- QKappa_attenuationext(i,j,ispec) = 10.d0
- Qmu_attenuationext(i,j,ispec) = 10.d0
+ QKappa_attenuationext(i,j,ispec) = 9999.d0
+ Qmu_attenuationext(i,j,ispec) = 9999.d0
else if((vsext(i,j,ispec) < TINYVAL) .and. (gravityext(i,j,ispec) < TINYVAL)) then
elastic(ispec) = .false.
poroelastic(ispec) = .false.
@@ -198,6 +207,17 @@
any_elastic = .true.
endif
+! attenuation is not implemented in acoustic (i.e. fluid) media for now, only in viscoelastic (i.e. solid) media
+ if(acoustic(ispec)) cycle
+
+! check that attenuation values entered by the user make sense
+ if((QKappa_attenuationext(i,j,ispec) <= 9998.999d0 .and. Qmu_attenuationext(i,j,ispec) > 9998.999d0) .or. &
+ (QKappa_attenuationext(i,j,ispec) > 9998.999d0 .and. Qmu_attenuationext(i,j,ispec) <= 9998.999d0)) stop &
+ 'need to have Qkappa and Qmu both above or both below 9999 for a given material; trick: use 9998 if you want to turn off one'
+
+! if no attenuation in that elastic element
+ if(QKappa_attenuationext(i,j,ispec) > 9998.999d0) cycle
+
call attenuation_model(N_SLS,QKappa_attenuationext(i,j,ispec),Qmu_attenuationext(i,j,ispec), &
f0_attenuation,tau_epsilon_nu1,tau_epsilon_nu2, &
inv_tau_sigma_nu1_sent,phi_nu1_sent,inv_tau_sigma_nu2_sent,phi_nu2_sent,Mu_nu1_sent,Mu_nu2_sent)
diff --git a/src/specfem2D/specfem2D.F90 b/src/specfem2D/specfem2D.F90
index d69f6ea..e2117d9 100644
--- a/src/specfem2D/specfem2D.F90
+++ b/src/specfem2D/specfem2D.F90
@@ -1372,7 +1372,7 @@
#ifdef USE_MPI
if(myrank == 0)then
if(time_stepping_scheme == 3) then
- stop 'mpi support for standard Runge Kutta scheme is not implemented'
+ stop 'MPI support for standard Runge-Kutta scheme is not implemented'
endif
endif
#endif
@@ -1423,13 +1423,35 @@
allocate(Mu_nu1(NGLLX,NGLLZ,nspec))
allocate(Mu_nu2(NGLLX,NGLLZ,nspec))
+! initialize to dummy values
+! convention to indicate that Q = 9999 in that element i.e. that there is no viscoelasticity in that element
+ inv_tau_sigma_nu1(:,:,:,:) = -1._CUSTOM_REAL
+ phi_nu1(:,:,:,:) = -1._CUSTOM_REAL
+ inv_tau_sigma_nu2(:,:,:,:) = -1._CUSTOM_REAL
+ phi_nu2(:,:,:,:) = -1._CUSTOM_REAL
+ Mu_nu1(:,:,:) = -1._CUSTOM_REAL
+ Mu_nu2(:,:,:) = -1._CUSTOM_REAL
+
! define the attenuation quality factors.
! they can be different for each element.
!! DK DK if needed in the future, here the quality factor could be different for each point
do ispec = 1,nspec
+
+! attenuation is not implemented in acoustic (i.e. fluid) media for now, only in viscoelastic (i.e. solid) media
+ if(acoustic(ispec)) cycle
+
+! check that attenuation values entered by the user make sense
+ if((QKappa_attenuation(kmato(ispec)) <= 9998.999d0 .and. Qmu_attenuation(kmato(ispec)) > 9998.999d0) .or. &
+ (QKappa_attenuation(kmato(ispec)) > 9998.999d0 .and. Qmu_attenuation(kmato(ispec)) <= 9998.999d0)) stop &
+ 'need to have Qkappa and Qmu both above or both below 9999 for a given material; trick: use 9998 if you want to turn off one'
+
+! if no attenuation in that elastic element
+ if(QKappa_attenuation(kmato(ispec)) > 9998.999d0) cycle
+
call attenuation_model(N_SLS,QKappa_attenuation(kmato(ispec)),Qmu_attenuation(kmato(ispec)), &
f0_attenuation,tau_epsilon_nu1,tau_epsilon_nu2,inv_tau_sigma_nu1_sent,phi_nu1_sent, &
inv_tau_sigma_nu2_sent,phi_nu2_sent,Mu_nu1_sent,Mu_nu2_sent)
+
do j = 1,NGLLZ
do i = 1,NGLLX
inv_tau_sigma_nu1(i,j,ispec,:) = inv_tau_sigma_nu1_sent(:)
@@ -1458,7 +1480,6 @@
poroelastcoef(3,1,n) = lambda + TWO*mu
already_shifted_velocity(n) = .true.
endif
-
endif
enddo
@@ -6303,7 +6324,7 @@ if(coupled_elastic_poro) then
rmemory_dux_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dx_prime,rmemory_duz_dz_prime, &
rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
- PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.false.,STACEY_BOUNDARY_CONDITIONS)
+ PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.false.,STACEY_BOUNDARY_CONDITIONS,acoustic)
if(SIMULATION_TYPE == 3)then
if(PML_BOUNDARY_CONDITIONS)then
@@ -6361,7 +6382,7 @@ if(coupled_elastic_poro) then
rmemory_dux_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dx_prime,rmemory_duz_dz_prime, &
rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
- .false.,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.true.,STACEY_BOUNDARY_CONDITIONS)
+ .false.,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.true.,STACEY_BOUNDARY_CONDITIONS,acoustic)
if(PML_BOUNDARY_CONDITIONS)then
do ispec = 1,nspec
More information about the CIG-COMMITS
mailing list