[cig-commits] r8566 - in seismo/2D/SPECFEM2D/trunk: .
UTILS/attenuation
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:57:01 PST 2007
Author: walter
Date: 2007-12-07 15:57:01 -0800 (Fri, 07 Dec 2007)
New Revision: 8566
Added:
seismo/2D/SPECFEM2D/trunk/attenuation_model.f90
Modified:
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/UTILS/attenuation/attenuation_code_2D.c
seismo/2D/SPECFEM2D/trunk/compute_energy.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
seismo/2D/SPECFEM2D/trunk/constants.h
seismo/2D/SPECFEM2D/trunk/constants_unstruct.h
seismo/2D/SPECFEM2D/trunk/construct_acoustic_surface.f90
seismo/2D/SPECFEM2D/trunk/create_color_image.f90
seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90
seismo/2D/SPECFEM2D/trunk/gll_library.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
modified code to be able to use any number (called N_SLS in constants.h) of
standard linear solids for attenuation (this was set to 2 before and
could not be changed). Also removed useless white spaces in some routines.
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-07 23:57:01 UTC (rev 8566)
@@ -2,7 +2,7 @@
# Makefile for SPECFEM2D version 5.2
#
# Dimitri Komatitsch, University of Pau, April 2007
-#
+#
SHELL=/bin/sh
O = obj
@@ -37,7 +37,7 @@
$O/lagrange_poly.o $O/gmat01.o $O/gll_library.o $O/plotgll.o $O/define_derivation_matrices.o\
$O/plotpost.o $O/locate_receivers.o $O/locate_source_force.o $O/compute_gradient_attenuation.o\
$O/specfem2D.o $O/write_seismograms.o $O/define_external_model.o $O/createnum_fast.o $O/createnum_slow.o\
- $O/define_shape_functions.o $O/create_color_image.o $O/compute_vector_field.o $O/compute_pressure.o\
+ $O/define_shape_functions.o $O/attenuation_model.o $O/create_color_image.o $O/compute_vector_field.o $O/compute_pressure.o\
$O/recompute_jacobian.o $O/compute_arrays_source.o $O/locate_source_moment_tensor.o $O/numerical_recipes.o\
$O/construct_acoustic_surface.o $O/assemble_MPI.o $O/compute_energy.o
@@ -112,6 +112,9 @@
$O/define_shape_functions.o: define_shape_functions.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/define_shape_functions.o define_shape_functions.f90
+$O/attenuation_model.o: attenuation_model.f90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/attenuation_model.o attenuation_model.f90
+
### use optimized compilation option for solver only
$O/specfem2D.o: specfem2D.F90 constants.h
${F90} $(FLAGS_NOCHECK) -c -o $O/specfem2D.o specfem2D.F90
Modified: seismo/2D/SPECFEM2D/trunk/UTILS/attenuation/attenuation_code_2D.c
===================================================================
--- seismo/2D/SPECFEM2D/trunk/UTILS/attenuation/attenuation_code_2D.c 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/UTILS/attenuation/attenuation_code_2D.c 2007-12-07 23:57:01 UTC (rev 8566)
@@ -19,8 +19,8 @@
main (argc,argv)
int argc; char **argv;
{
- int xmgr, n, i, j, plot;
- float Q_kappa, Q_s, Q_kappa_m, Q_s_m;
+ int xmgr, n, i, j, plot, nu;
+ double Q_s, target_Qp, target_Qs;
double f1, f2, Q, om0, Omega;
double a, b;
double kappa, mu, kappa0, mu0, kappaR, muR;
@@ -29,10 +29,14 @@
void constant_Q2_sub(),plot_modulus();
void free_dvector();
- printf("value of Q_p or Q_s: ");
- scanf("%f",&Q_s);
- printf("%f\n",Q_s);
+ printf("target value of Qp: ");
+ scanf("%lf",&target_Qp);
+ printf("%lf\n",target_Qp);
+ printf("target value of Qs: ");
+ scanf("%lf",&target_Qs);
+ printf("%lf\n",target_Qs);
+
printf("shortest frequency (Hz): ");
scanf("%lf",&f1);
printf("%lf\n",f1);
@@ -51,65 +55,67 @@
if (f2 < f1) {
printf("T2 > T1\n");
- exit;
- }
- if (Q < 0.0) {
- printf("Q < 0\n");
- exit;
- }
+ exit; }
+
+ if (target_Qp <= 0.0001) {
+ printf("Qp cannot be negative or null\n");
+ exit; }
+
+ if (target_Qs <= 0.0001) {
+ printf("Qs cannot be negative or null\n");
+ exit; }
+
if (n < 1) {
printf("n < 1\n");
- exit;
- }
+ exit; }
- tau_s = dvector(1, n);
- tau_e = dvector(1, n);
-
- Q_kappa_m = Q_s_m = 0.0;
om0 = PI2 * pow(10.0, 0.5 * (log10(f1) + log10(f2)));
-/* DK DK printf("\n\n! central frequency: %25.15f mHz\n\n", 1.0E+03 * om0 / PI2); */
- plot=0;
+ printf("\n! put this in file constants.h\n\n");
-/* DK DK removed for Qmu only in the Earth
- if (Q_kappa != Q_kappa_m) {
- printf("\ntarget Q_kappa: %6.2f\n\n", Q_kappa);
- constant_Q2_sub(f1, f2, n, (double) Q_kappa, tau_s, tau_e, xmgr);
- Q_kappa_m = Q_kappa;
- a = 1.0;
- b = 0.0;
- for (i = 1; i <= n; i++) {
- a -= om0 * om0 * tau_e[i] * (tau_e[i] - tau_s[i]) /
- (1.0 + om0 * om0 * tau_e[i] * tau_e[i]);
- b += om0 * (tau_e[i] - tau_s[i]) /
- (1.0 + om0 * om0 * tau_e[i] * tau_e[i]);
- }
- }
-*/
+ printf("! number of standard linear solids for attenuation\n");
+ printf(" integer, parameter :: N_SLS = %d\n\n",n);
- printf("\n! frequency range: %f Hz - %f Hz\n\n", f1 , f2);
- printf("! central frequency in log scale in Hz = %20.15f\n",om0 / PI2);
+ printf("! put this in file attenuation_model.f90\n\n");
- constant_Q2_sub(f1, f2, n, (double) Q_s, tau_s, tau_e, xmgr);
+ printf("! frequency range: %lf Hz - %lf Hz\n", f1 , f2);
+ printf("! central frequency in log scale in Hz = %20.15f\n",om0 / PI2);
-/* DK DK converted to Fortran90 output */
+ printf("! target constant attenuation factor Qp = %20.10lf\n", target_Qp);
+ printf("! target constant attenuation factor Qs = %20.10lf\n\n", target_Qs);
- printf("! target constant attenuation factor Q = %20.10f\n\n", Q_s);
+ printf("! tau_sigma evenly spaced in log frequency, do not depend on value of Q\n\n");
- printf("! tau sigma evenly spaced in log frequency, do not depend on value of Q\n");
- for (i = 1; i <= n; i++) {
- printf(" double precision, parameter :: tau_sigma_mech%1d = %30.20fd0\n", i, tau_s[i]);
- }
- printf("\n");
+ plot = 0;
- for (i = 1; i <= n; i++) {
- printf(" double precision, parameter :: tau_epsilon_mech%1d = %30.20fd0\n", i, tau_e[i] );
- }
- printf("\n");
+/* loop on the Qp dilatation mode (nu = 1) and Qs shear mode (nu = 2) */
+ for (nu = 1; nu <= 2; nu++) {
- free_dvector(tau_s, 1, n);
- free_dvector(tau_e, 1, n);
+/* assign Qp or Qs to generic variable Q_s which is used for the calculations */
+ if (nu == 1) { Q_s = target_Qp ; }
+ if (nu == 2) { Q_s = target_Qs ; }
+ tau_s = dvector(1, n);
+ tau_e = dvector(1, n);
+
+ constant_Q2_sub(f1, f2, n, Q_s, tau_s, tau_e, xmgr);
+
+/* output in Fortran90 format */
+ for (i = 1; i <= n; i++) {
+ printf(" tau_sigma_nu%d(%1d) = %30.20lfd0\n", nu, i, tau_s[i]);
+ }
+ printf("\n");
+
+ for (i = 1; i <= n; i++) {
+ printf(" tau_epsilon_nu%d(%1d) = %30.20lfd0\n", nu, i, tau_e[i]);
+ }
+ printf("\n");
+
+ free_dvector(tau_s, 1, n);
+ free_dvector(tau_e, 1, n);
+
+ }
+
}
void plot_modulus(f1, f2, n, m, mR, Q, tau_e, tau_s ,xmgr)
@@ -905,7 +911,7 @@
double *x1, *x2;
double *gradient, **hessian;
double *dvector(), **dmatrix();
- void print_model(), derivatives();
+ void derivatives();
void initialize(), invert();
void free_dvector(), free_dmatrix();
@@ -942,20 +948,11 @@
free_dvector(gradient, 1, n);
free_dmatrix(hessian, 1, n, 1, n);
- print_model(f1, f2, n, Q, x1, x2, xmgr);
-
-/* DK DK printf("! frequency range: %f -- %f mHz\n", f1 * 1.0E+03, f2 * 1.0E+03);
- printf("! period range: %f -- %f s\n", 1./f2 , 1./f1 );
- printf("! desired Q: %f\n", Q);
- printf("! number of relaxation mechanisms: %d\n", n); */
for (i = 1; i <= n; i++) {
tau_e[i]=x1[i] + x2[i];
-/* DK DK printf("tau_e[%d]: %e\n", i, x1[i] + x2[i]); */
}
- printf("\n");
for (i = 1; i <= n; i++) {
tau_s[i]=x2[i];
-/* DK DK printf("tau_s[%d]: %e\n", i, x2[i]); */
}
free_dvector(x1, 1, n);
@@ -1002,18 +999,6 @@
x2[i] = tau_s[i];
}
-/* DK DK suppressed verbose output
-printf("initial stress and strain relaxation times: \n\n");
-for (i = 1; i <= n; i++) {
- printf("tau_e[%d]: %e\n", i, x1[i] + x2[i]);
-}
-printf("\n");
-for (i = 1; i <= n; i++) {
- printf("tau_s[%d]: %e\n", i, x2[i]);
-}
-printf("\n"); */
-
-
free_dvector(tau_e, 1, n);
free_dvector(tau_s, 1, n);
}
@@ -1053,64 +1038,7 @@
return pnlt;
}
-void print_model(f1, f2, n, Q, x1, x2, xmgr)
- int n, xmgr;
- double f1, f2, Q, *x1, *x2;
-{
-int pid, i;
-double exp1, exp2, dexp, expo;
-double f, omega;
-double tau_e, tau_s, a, b, Q_omega;
-char strng[180];
-int getpid(), system();
-FILE *fp_q, *fp_q_approx;
-pid = getpid();
-sprintf(strng, "q%1d", pid);
-if((fp_q=fopen(strng,"w"))==NULL) {
- puts("cannot open file\n");
- exit;
-}
-sprintf(strng, "q_approx%1d", pid);
-if((fp_q_approx=fopen(strng,"w"))==NULL) {
- puts("cannot open file\n");
- exit;
-}
-
-exp1 = log10(f1) - 2.0;
-exp2 = log10(f2) + 2.0;
-dexp = (exp2 - exp1) / 100.0;
-for (expo = exp1; expo <= exp2; expo += dexp) {
- f = pow(10.0, expo);
- omega = PI2 * f;
- a = (double) (1 - n);
- b = 0.0;
- for (i = 1; i <= n; i++) {
- tau_e = x1[i] + x2[i];
- tau_s = x2[i];
- a += (1.0 + omega * omega * tau_e * tau_s) /
- (1.0 + omega * omega * tau_s * tau_s);
- b += omega * (tau_e - tau_s) /
- (1.0 + omega * omega * tau_s * tau_s);
- }
- Q_omega = a / b;
- if (omega >= PI2 * f1 && omega <= PI2 * f2) {
- fprintf(fp_q, "%f %f\n", f, Q);
- fprintf(fp_q_approx, "%f %f\n", f, Q_omega);
- }
-}
-fclose(fp_q);
-fclose(fp_q_approx);
-
-/* DK DK added option to avoid calling Xmgr */
-if(xmgr == 1) {
- sprintf(strng, "xmgr q%1d q_approx%1d", pid, pid);
- system(strng);
- sprintf(strng, "rm q%1d q_approx%1d", pid, pid, pid);
- system(strng);
-}
-}
-
void derivatives(f1, f2, n, Q, x1, x2, gradient, hessian)
int n;
double f1, f2, Q, *x1, *x2;
Added: seismo/2D/SPECFEM2D/trunk/attenuation_model.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/attenuation_model.f90 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/attenuation_model.f90 2007-12-07 23:57:01 UTC (rev 8566)
@@ -0,0 +1,114 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.2
+! ------------------------------
+!
+! Dimitri Komatitsch
+! University of Pau, France
+!
+! (c) April 2007
+!
+!========================================================================
+
+ subroutine attenuation_model(inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2)
+
+! define the attenuation constants
+
+ implicit none
+
+ include "constants.h"
+
+ double precision, dimension(N_SLS) :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
+ double precision :: Mu_nu1,Mu_nu2
+
+ integer :: i_sls
+
+ double precision, dimension(N_SLS) :: tau_epsilon_nu1,tau_sigma_nu1,tau_epsilon_nu2,tau_sigma_nu2
+
+! attenuation constants for standard linear solids
+
+! nu1 is the dilatation mode
+! nu2 is the shear mode
+
+! array index (1) is the first standard linear solid, (2) is the second etc.
+
+! from J. M. Carcione, Seismic modeling in viscoelastic media, Geophysics,
+! vol. 58(1), p. 110-120 (1993) for two memory-variable mechanisms (page 112).
+! Beware: these values implement specific values of the quality factors:
+! Qp approximately equal to 13 and Qs approximately equal to 10,
+! which means very high attenuation, see that paper for details.
+! tau_epsilon_nu1(1) = 0.0334d0
+! tau_sigma_nu1(1) = 0.0303d0
+! tau_epsilon_nu2(1) = 0.0352d0
+! tau_sigma_nu2(1) = 0.0287d0
+
+! tau_epsilon_nu1(2) = 0.0028d0
+! tau_sigma_nu1(2) = 0.0025d0
+! tau_epsilon_nu2(2) = 0.0029d0
+! tau_sigma_nu2(2) = 0.0024d0
+
+! from J. M. Carcione, D. Kosloff and R. Kosloff, Wave propagation simulation
+! in a linear viscoelastic medium, Geophysical Journal International,
+! vol. 95, p. 597-611 (1988) for two memory-variable mechanisms (page 604).
+! Beware: these values implement specific values of the quality factors:
+! Qp approximately equal to 27 and Qs approximately equal to 20,
+! which means very high attenuation, see that paper for details.
+ tau_epsilon_nu1(1) = 0.0325305d0
+ tau_sigma_nu1(1) = 0.0311465d0
+ tau_epsilon_nu2(1) = 0.0332577d0
+ tau_sigma_nu2(1) = 0.0304655d0
+
+ tau_epsilon_nu1(2) = 0.0032530d0
+ tau_sigma_nu1(2) = 0.0031146d0
+ tau_epsilon_nu2(2) = 0.0033257d0
+ tau_sigma_nu2(2) = 0.0030465d0
+
+! values for Paul Cristini for fluid-solid ocean acoustics simulations
+
+! frequency range: 3.000000 Hz - 20.000000 Hz
+! central frequency in log scale in Hz = 7.745966692414834
+! target constant attenuation factor Q = 136.4376068115
+
+! tau sigma evenly spaced in log frequency, do not depend on value of Q
+! tau_sigma_nu1(1) = 0.05305164769729849711d0
+! tau_sigma_nu1(2) = 0.00795774715459477387d0
+
+! tau_epsilon_nu1(1) = 0.05361741010508015715d0
+! tau_epsilon_nu1(2) = 0.00804740719550106794d0
+
+! frequency range: 1.500000 Hz - 18.000000 Hz
+! central frequency in log scale in Hz = 5.196152422706633
+! target constant attenuation factor Q = 136.4376068115
+
+! tau sigma evenly spaced in log frequency, do not depend on value of Q
+! tau_sigma_nu1(1) = 0.10610329539459699422d0
+! tau_sigma_nu1(2) = 0.00884194128288308401d0
+
+! tau_epsilon_nu1(1) = 0.10754721280605997191d0
+! tau_epsilon_nu1(2) = 0.00895488050110176612d0
+
+! tau_epsilon_nu2(1) = tau_epsilon_nu1(1)
+! tau_epsilon_nu2(2) = tau_epsilon_nu1(2)
+! tau_sigma_nu2(1) = tau_sigma_nu1(1)
+! tau_sigma_nu2(2) = tau_sigma_nu1(2)
+
+!
+!--- other constants computed from the parameters above, do not modify
+!
+ inv_tau_sigma_nu1(:) = ONE / tau_sigma_nu1(:)
+ inv_tau_sigma_nu2(:) = ONE / tau_sigma_nu2(:)
+
+ phi_nu1(:) = (ONE - tau_epsilon_nu1(:)/tau_sigma_nu1(:)) / tau_sigma_nu1(:)
+ phi_nu2(:) = (ONE - tau_epsilon_nu2(:)/tau_sigma_nu2(:)) / tau_sigma_nu2(:)
+
+ Mu_nu1 = ONE
+ Mu_nu2 = ONE
+
+ do i_sls = 1,N_SLS
+ Mu_nu1 = Mu_nu1 - (ONE - tau_epsilon_nu1(i_sls)/tau_sigma_nu1(i_sls))
+ Mu_nu2 = Mu_nu2 - (ONE - tau_epsilon_nu2(i_sls)/tau_sigma_nu2(i_sls))
+ enddo
+
+ end subroutine attenuation_model
+
Modified: seismo/2D/SPECFEM2D/trunk/compute_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2007-12-07 23:57:01 UTC (rev 8566)
@@ -15,8 +15,8 @@
xix,xiz,gammax,gammaz,jacobian,ibool,elastic,hprime_xx,hprime_zz, &
nspec,npoin,assign_external_model,it,deltat,t0,kmato,elastcoef,density, &
vpext,vsext,rhoext,wxgll,wzgll,numat, &
- pressure_element,vector_field_element,e1_mech1,e11_mech1,e1_mech2,e11_mech2, &
- potential_dot_acoustic,potential_dot_dot_acoustic,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON)
+ pressure_element,vector_field_element,e1,e11, &
+ potential_dot_acoustic,potential_dot_dot_acoustic,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
! compute kinetic and potential energy in the solid (acoustic elements are excluded)
@@ -30,7 +30,8 @@
! pressure in an element
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: pressure_element
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: e1_mech1,e11_mech1,e1_mech2,e11_mech2
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11
+ double precision :: Mu_nu1,Mu_nu2
real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_acoustic,potential_dot_dot_acoustic
@@ -170,8 +171,8 @@
! compute pressure in this element
call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,elastic, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1_mech1,e11_mech1, &
- e1_mech2,e11_mech2,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON)
+ numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
+ TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
! compute velocity vector field in this element
call compute_vector_one_element(vector_field_element,potential_dot_acoustic,veloc_elastic,elastic, &
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2007-12-07 23:57:01 UTC (rev 8566)
@@ -16,10 +16,10 @@
initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,angleforce,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,elastic,codeabs, &
accel_elastic,veloc_elastic,displ_elastic,density,elastcoef,xix,xiz,gammax,gammaz, &
- jacobian,vpext,vsext,rhoext,source_time_function,sourcearray,e1_mech1,e11_mech1, &
- e13_mech1,e1_mech2,e11_mech2,e13_mech2,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
+ jacobian,vpext,vsext,rhoext,source_time_function,sourcearray,e1,e11, &
+ e13,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll)
+ hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2)
! compute forces for the elastic elements
@@ -48,8 +48,13 @@
real(kind=CUSTOM_REAL), dimension(NSTEP) :: source_time_function
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLZ) :: sourcearray
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11,e13
+ double precision, dimension(N_SLS) :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
+ double precision :: Mu_nu1,Mu_nu2
+ real(kind=CUSTOM_REAL) :: e1_sum,e11_sum,e13_sum
+ integer :: i_sls
+
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: &
- e1_mech1,e11_mech1,e13_mech1,e1_mech2,e11_mech2,e13_mech2, &
dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1
! derivatives of Lagrange polynomials
@@ -165,12 +170,20 @@
! add the memory variables using the relaxed parameters (Carcione 1993, page 111)
! beware: there is a bug in Carcione's equation (2c) for sigma_zz, we fixed it in the code below
- sigma_xx = sigma_xx + (lambdal_relaxed + mul_relaxed)* &
- (e1_mech1(i,j,ispec) + e1_mech2(i,j,ispec)) + TWO * mul_relaxed * (e11_mech1(i,j,ispec) + e11_mech2(i,j,ispec))
- sigma_xz = sigma_xz + mul_relaxed * (e13_mech1(i,j,ispec) + e13_mech2(i,j,ispec))
- sigma_zz = sigma_zz + (lambdal_relaxed + mul_relaxed)* &
- (e1_mech1(i,j,ispec) + e1_mech2(i,j,ispec)) - TWO * mul_relaxed * (e11_mech1(i,j,ispec) + e11_mech2(i,j,ispec))
+ e1_sum = 0._CUSTOM_REAL
+ e11_sum = 0._CUSTOM_REAL
+ e13_sum = 0._CUSTOM_REAL
+ do i_sls = 1,N_SLS
+ e1_sum = e1_sum + e1(i,j,ispec,i_sls)
+ e11_sum = e11_sum + e11(i,j,ispec,i_sls)
+ e13_sum = e13_sum + e13(i,j,ispec,i_sls)
+ enddo
+
+ sigma_xx = sigma_xx + (lambdal_relaxed + mul_relaxed) * e1_sum + TWO * mul_relaxed * e11_sum
+ sigma_xz = sigma_xz + mul_relaxed * e13_sum
+ sigma_zz = sigma_zz + (lambdal_relaxed + mul_relaxed) * e1_sum - TWO * mul_relaxed * e11_sum
+
else
! no attenuation
@@ -490,90 +503,53 @@
theta_n = dux_dxl_n(i,j,ispec) + duz_dzl_n(i,j,ispec)
theta_np1 = dux_dxl_np1(i,j,ispec) + duz_dzl_np1(i,j,ispec)
-! evolution e1_mech1
- Un = e1_mech1(i,j,ispec)
- tauinv = - inv_tau_sigma_nu1_mech1
- tauinvsquare = tauinv * tauinv
- tauinvcube = tauinvsquare * tauinv
- tauinvUn = tauinv * Un
- Sn = theta_n * phi_nu1_mech1
- Snp1 = theta_np1 * phi_nu1_mech1
- Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
- twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
- fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
- deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e1_mech1(i,j,ispec) = Unp1
+! loop on all the standard linear solids
+ do i_sls = 1,N_SLS
-! evolution e1_mech2
- Un = e1_mech2(i,j,ispec)
- tauinv = - inv_tau_sigma_nu1_mech2
+! evolution e1
+ Un = e1(i,j,ispec,i_sls)
+ tauinv = - inv_tau_sigma_nu1(i_sls)
tauinvsquare = tauinv * tauinv
tauinvcube = tauinvsquare * tauinv
tauinvUn = tauinv * Un
- Sn = theta_n * phi_nu1_mech2
- Snp1 = theta_np1 * phi_nu1_mech2
+ Sn = theta_n * phi_nu1(i_sls)
+ Snp1 = theta_np1 * phi_nu1(i_sls)
Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e1_mech2(i,j,ispec) = Unp1
+ e1(i,j,ispec,i_sls) = Unp1
-! evolution e11_mech1
- Un = e11_mech1(i,j,ispec)
- tauinv = - inv_tau_sigma_nu2_mech1
+! evolution e11
+ Un = e11(i,j,ispec,i_sls)
+ tauinv = - inv_tau_sigma_nu2(i_sls)
tauinvsquare = tauinv * tauinv
tauinvcube = tauinvsquare * tauinv
tauinvUn = tauinv * Un
- Sn = (dux_dxl_n(i,j,ispec) - theta_n/TWO) * phi_nu2_mech1
- Snp1 = (dux_dxl_np1(i,j,ispec) - theta_np1/TWO) * phi_nu2_mech1
+ Sn = (dux_dxl_n(i,j,ispec) - theta_n/TWO) * phi_nu2(i_sls)
+ Snp1 = (dux_dxl_np1(i,j,ispec) - theta_np1/TWO) * phi_nu2(i_sls)
Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e11_mech1(i,j,ispec) = Unp1
+ e11(i,j,ispec,i_sls) = Unp1
-! evolution e11_mech2
- Un = e11_mech2(i,j,ispec)
- tauinv = - inv_tau_sigma_nu2_mech2
+! evolution e13
+ Un = e13(i,j,ispec,i_sls)
+ tauinv = - inv_tau_sigma_nu2(i_sls)
tauinvsquare = tauinv * tauinv
tauinvcube = tauinvsquare * tauinv
tauinvUn = tauinv * Un
- Sn = (dux_dxl_n(i,j,ispec) - theta_n/TWO) * phi_nu2_mech2
- Snp1 = (dux_dxl_np1(i,j,ispec) - theta_np1/TWO) * phi_nu2_mech2
+ Sn = (dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)) * phi_nu2(i_sls)
+ Snp1 = (dux_dzl_np1(i,j,ispec) + duz_dxl_np1(i,j,ispec)) * phi_nu2(i_sls)
Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e11_mech2(i,j,ispec) = Unp1
+ e13(i,j,ispec,i_sls) = Unp1
-! evolution e13_mech1
- Un = e13_mech1(i,j,ispec)
- tauinv = - inv_tau_sigma_nu2_mech1
- tauinvsquare = tauinv * tauinv
- tauinvcube = tauinvsquare * tauinv
- tauinvUn = tauinv * Un
- Sn = (dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)) * phi_nu2_mech1
- Snp1 = (dux_dzl_np1(i,j,ispec) + duz_dxl_np1(i,j,ispec)) * phi_nu2_mech1
- Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
- twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
- fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
- deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e13_mech1(i,j,ispec) = Unp1
+ enddo
-! evolution e13_mech2
- Un = e13_mech2(i,j,ispec)
- tauinv = - inv_tau_sigma_nu2_mech2
- tauinvsquare = tauinv * tauinv
- tauinvcube = tauinvsquare * tauinv
- tauinvUn = tauinv * Un
- Sn = (dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)) * phi_nu2_mech2
- Snp1 = (dux_dzl_np1(i,j,ispec) + duz_dxl_np1(i,j,ispec)) * phi_nu2_mech2
- Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
- twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
- fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
- deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e13_mech2(i,j,ispec) = Unp1
-
enddo
enddo
enddo
Modified: seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2007-12-07 23:57:01 UTC (rev 8566)
@@ -13,8 +13,8 @@
subroutine compute_pressure_whole_medium(potential_dot_dot_acoustic,displ_elastic,elastic,vector_field_display, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,elastcoef,vpext,vsext,rhoext,e1_mech1,e11_mech1, &
- e1_mech2,e11_mech2,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON)
+ numat,kmato,density,elastcoef,vpext,vsext,rhoext,e1,e11, &
+ TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
! compute pressure in acoustic elements and in elastic elements
@@ -44,7 +44,8 @@
logical :: assign_external_model,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: e1_mech1,e11_mech1,e1_mech2,e11_mech2
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11
+ double precision :: Mu_nu1,Mu_nu2
! local variables
integer :: i,j,ispec,iglob
@@ -58,8 +59,8 @@
! compute pressure in this element
call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,elastic, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1_mech1,e11_mech1, &
- e1_mech2,e11_mech2,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON)
+ numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
+ TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
! use vector_field_display as temporary storage, store pressure in its second component
do j = 1,NGLLZ
@@ -79,8 +80,8 @@
subroutine compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,elastic, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1_mech1,e11_mech1, &
- e1_mech2,e11_mech2,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON)
+ numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
+ TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
! compute pressure in acoustic elements and in elastic elements
@@ -112,7 +113,10 @@
logical :: assign_external_model,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: e1_mech1,e11_mech1,e1_mech2,e11_mech2
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11
+ real(kind=CUSTOM_REAL) :: e1_sum,e11_sum
+ double precision :: Mu_nu1,Mu_nu2
+ integer :: i_sls
! local variables
integer :: i,j,k,iglob
@@ -217,11 +221,17 @@
! add the memory variables using the relaxed parameters (Carcione 1993, page 111)
! beware: there is a bug in Carcione's equation (2c) for sigma_zz, we fixed it in the code below
- sigma_xx = sigma_xx + (lambdal_relaxed + mul_relaxed)* &
- (e1_mech1(i,j,ispec) + e1_mech2(i,j,ispec)) + TWO * mul_relaxed * (e11_mech1(i,j,ispec) + e11_mech2(i,j,ispec))
- sigma_zz = sigma_zz + (lambdal_relaxed + mul_relaxed)* &
- (e1_mech1(i,j,ispec) + e1_mech2(i,j,ispec)) - TWO * mul_relaxed * (e11_mech1(i,j,ispec) + e11_mech2(i,j,ispec))
+ e1_sum = 0._CUSTOM_REAL
+ e11_sum = 0._CUSTOM_REAL
+ do i_sls = 1,N_SLS
+ e1_sum = e1_sum + e1(i,j,ispec,i_sls)
+ e11_sum = e11_sum + e11(i,j,ispec,i_sls)
+ enddo
+
+ sigma_xx = sigma_xx + (lambdal_relaxed + mul_relaxed) * e1_sum + TWO * mul_relaxed * e11_sum
+ sigma_zz = sigma_zz + (lambdal_relaxed + mul_relaxed) * e1_sum - TWO * mul_relaxed * e11_sum
+
else
! no attenuation
Modified: seismo/2D/SPECFEM2D/trunk/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants.h 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/constants.h 2007-12-07 23:57:01 UTC (rev 8566)
@@ -14,6 +14,9 @@
! the code does NOT work if NGLLZ /= NGLLX because it then cannot handle a non-structured mesh
integer, parameter :: NGLLZ = NGLLX
+! number of standard linear solids for attenuation
+ integer, parameter :: N_SLS = 2
+
! compute and output acoustic and elastic energy (slows down the code significantly)
logical, parameter :: OUTPUT_ENERGY = .false.
@@ -100,93 +103,7 @@
!-----------------------------------------------------------------------
-! attenuation constants for standard linear solids
-! nu1 is the dilatation mode
-! nu2 is the shear mode
-! mech1 is the first standard linear solid, mech2 is the second
-
-! from J. M. Carcione, Seismic modeling in viscoelastic media, Geophysics,
-! vol. 58(1), p. 110-120 (1993) for two memory-variable mechanisms (page 112).
-! Beware: these values implement specific values of the quality factors:
-! Qp approximately equal to 13 and Qs approximately equal to 10,
-! which means very high attenuation, see that paper for details.
-! double precision, parameter :: tau_epsilon_nu1_mech1 = 0.0334d0
-! double precision, parameter :: tau_sigma_nu1_mech1 = 0.0303d0
-! double precision, parameter :: tau_epsilon_nu2_mech1 = 0.0352d0
-! double precision, parameter :: tau_sigma_nu2_mech1 = 0.0287d0
-
-! double precision, parameter :: tau_epsilon_nu1_mech2 = 0.0028d0
-! double precision, parameter :: tau_sigma_nu1_mech2 = 0.0025d0
-! double precision, parameter :: tau_epsilon_nu2_mech2 = 0.0029d0
-! double precision, parameter :: tau_sigma_nu2_mech2 = 0.0024d0
-
-! from J. M. Carcione, D. Kosloff and R. Kosloff, Wave propagation simulation
-! in a linear viscoelastic medium, Geophysical Journal International,
-! vol. 95, p. 597-611 (1988) for two memory-variable mechanisms (page 604).
-! Beware: these values implement specific values of the quality factors:
-! Qp approximately equal to 27 and Qs approximately equal to 20,
-! which means very high attenuation, see that paper for details.
- double precision, parameter :: tau_epsilon_nu1_mech1 = 0.0325305d0
- double precision, parameter :: tau_sigma_nu1_mech1 = 0.0311465d0
- double precision, parameter :: tau_epsilon_nu2_mech1 = 0.0332577d0
- double precision, parameter :: tau_sigma_nu2_mech1 = 0.0304655d0
-
- double precision, parameter :: tau_epsilon_nu1_mech2 = 0.0032530d0
- double precision, parameter :: tau_sigma_nu1_mech2 = 0.0031146d0
- double precision, parameter :: tau_epsilon_nu2_mech2 = 0.0033257d0
- double precision, parameter :: tau_sigma_nu2_mech2 = 0.0030465d0
-
-! values for Paul Cristini for fluid-solid ocean acoustics simulations
-
-! frequency range: 3.000000 Hz - 20.000000 Hz
-! central frequency in log scale in Hz = 7.745966692414834
-! target constant attenuation factor Q = 136.4376068115
-
-! tau sigma evenly spaced in log frequency, do not depend on value of Q
-! double precision, parameter :: tau_sigma_nu1_mech1 = 0.05305164769729849711d0
-! double precision, parameter :: tau_sigma_nu1_mech2 = 0.00795774715459477387d0
-
-! double precision, parameter :: tau_epsilon_nu1_mech1 = 0.05361741010508015715d0
-! double precision, parameter :: tau_epsilon_nu1_mech2 = 0.00804740719550106794d0
-
-! frequency range: 1.500000 Hz - 18.000000 Hz
-! central frequency in log scale in Hz = 5.196152422706633
-! target constant attenuation factor Q = 136.4376068115
-
-! tau sigma evenly spaced in log frequency, do not depend on value of Q
-! double precision, parameter :: tau_sigma_nu1_mech1 = 0.10610329539459699422d0
-! double precision, parameter :: tau_sigma_nu1_mech2 = 0.00884194128288308401d0
-
-! double precision, parameter :: tau_epsilon_nu1_mech1 = 0.10754721280605997191d0
-! double precision, parameter :: tau_epsilon_nu1_mech2 = 0.00895488050110176612d0
-
-! double precision, parameter :: tau_epsilon_nu2_mech1 = tau_epsilon_nu1_mech1
-! double precision, parameter :: tau_epsilon_nu2_mech2 = tau_epsilon_nu1_mech2
-! double precision, parameter :: tau_sigma_nu2_mech1 = tau_sigma_nu1_mech1
-! double precision, parameter :: tau_sigma_nu2_mech2 = tau_sigma_nu1_mech2
-
!
-!--- other constants computed from the parameters above, do not modify
-!
-
- double precision, parameter :: inv_tau_sigma_nu1_mech1 = ONE / tau_sigma_nu1_mech1
- double precision, parameter :: inv_tau_sigma_nu2_mech1 = ONE / tau_sigma_nu2_mech1
- double precision, parameter :: inv_tau_sigma_nu1_mech2 = ONE / tau_sigma_nu1_mech2
- double precision, parameter :: inv_tau_sigma_nu2_mech2 = ONE / tau_sigma_nu2_mech2
-
- double precision, parameter :: phi_nu1_mech1 = (ONE - tau_epsilon_nu1_mech1/tau_sigma_nu1_mech1) / tau_sigma_nu1_mech1
- double precision, parameter :: phi_nu2_mech1 = (ONE - tau_epsilon_nu2_mech1/tau_sigma_nu2_mech1) / tau_sigma_nu2_mech1
- double precision, parameter :: phi_nu1_mech2 = (ONE - tau_epsilon_nu1_mech2/tau_sigma_nu1_mech2) / tau_sigma_nu1_mech2
- double precision, parameter :: phi_nu2_mech2 = (ONE - tau_epsilon_nu2_mech2/tau_sigma_nu2_mech2) / tau_sigma_nu2_mech2
-
- double precision, parameter :: Mu_nu1 = ONE - (ONE - tau_epsilon_nu1_mech1/tau_sigma_nu1_mech1) &
- - (ONE - tau_epsilon_nu1_mech2/tau_sigma_nu1_mech2)
- double precision, parameter :: Mu_nu2 = ONE - (ONE - tau_epsilon_nu2_mech1/tau_sigma_nu2_mech1) &
- - (ONE - tau_epsilon_nu2_mech2/tau_sigma_nu2_mech2)
-
-!-----------------------------------------------------------------------
-
-!
! anisotropic copper crystal (cubic symmetry)
!
Modified: seismo/2D/SPECFEM2D/trunk/constants_unstruct.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants_unstruct.h 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/constants_unstruct.h 2007-12-07 23:57:01 UTC (rev 8566)
@@ -1,5 +1,5 @@
-! Number of nodes per elements.
+! Number of nodes per elements.
integer, parameter :: ESIZE = 4
! Max number of neighbours per elements.
Modified: seismo/2D/SPECFEM2D/trunk/construct_acoustic_surface.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/construct_acoustic_surface.f90 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/construct_acoustic_surface.f90 2007-12-07 23:57:01 UTC (rev 8566)
@@ -1,5 +1,5 @@
subroutine construct_acoustic_surface ( nspec, ngnod, knods, nsurface, surface, tab_surface )
-
+
implicit none
integer, intent(in) :: nspec
@@ -8,35 +8,35 @@
integer, intent(in) :: nsurface
integer, dimension(4,nsurface), intent(in) :: surface
integer, dimension(5,nsurface), intent(out) :: tab_surface
-
+
integer :: i, k
integer :: ixmin, ixmax
integer :: izmin, izmax
integer, dimension(ngnod) :: n
integer :: e1, e2
integer :: type
-
-
+
+
do i = 1, nsurface
tab_surface(1,i) = surface(1,i)
type = surface(2,i)
e1 = surface(3,i)
e2 = surface(4,i)
do k = 1, ngnod
- n(k) = knods(k,tab_surface(1,i))
+ n(k) = knods(k,tab_surface(1,i))
end do
-
+
call get_acoustic_edge ( ngnod, n, type, e1, e2, ixmin, ixmax, izmin, izmax )
-
+
tab_surface(2,i) = ixmin
tab_surface(3,i) = ixmax
tab_surface(4,i) = izmin
tab_surface(5,i) = izmax
-
+
end do
-
+
end subroutine construct_acoustic_surface
@@ -51,8 +51,8 @@
integer, dimension(ngnod), intent(in) :: n
integer, intent(in) :: type, e1, e2
integer, intent(out) :: ixmin, ixmax, izmin, izmax
-
-
+
+
if ( type == 1 ) then
if ( e1 == n(1) ) then
ixmin = 1
@@ -78,20 +78,20 @@
izmin = NGLLZ
izmax = NGLLZ
end if
-
+
else
if ( e1 == n(1) ) then
ixmin = 1
izmin = 1
- if ( e2 == n(2) ) then
+ if ( e2 == n(2) ) then
ixmax = NGLLX
izmax = 1
-
+
end if
if ( e2 == n(4) ) then
ixmax = 1
izmax = NGLLZ
-
+
end if
end if
if ( e1 == n(2) ) then
@@ -100,13 +100,13 @@
if ( e2 == n(3) ) then
ixmax = NGLLX
izmax = NGLLZ
-
+
end if
if ( e2 == n(1) ) then
ixmax = ixmin
ixmin = 1
izmax = 1
-
+
end if
end if
if ( e1 == n(3) ) then
@@ -122,7 +122,7 @@
ixmax = NGLLX
izmax = izmin
izmin = 1
-
+
end if
end if
if ( e1 == n(4) ) then
@@ -132,17 +132,17 @@
ixmax = 1
izmax = izmin
izmin = 1
-
+
end if
if ( e2 == n(3) ) then
ixmax = NGLLX
izmax = NGLLZ
-
+
end if
end if
end if
-
+
end subroutine get_acoustic_edge
Modified: seismo/2D/SPECFEM2D/trunk/create_color_image.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/create_color_image.f90 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/create_color_image.f90 2007-12-07 23:57:01 UTC (rev 8566)
@@ -215,7 +215,7 @@
! open image file and create system command to convert image to more convenient format
write(system_command,"('cd OUTPUT_FILES ; convert image',i7.7,'.pnm image',i7.7,'.gif')") it,it
-
+
! call the system to convert image to GIF
call system(system_command)
Modified: seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90 2007-12-07 23:57:01 UTC (rev 8566)
@@ -36,9 +36,9 @@
!---
integer :: ispec_acoustic_surface,ispec,i,j,iglob
-
+
do ispec_acoustic_surface = 1, nelem_acoustic_surface
-
+
ispec = acoustic_surface(1,ispec_acoustic_surface)
do j = acoustic_surface(4,ispec_acoustic_surface), acoustic_surface(5,ispec_acoustic_surface)
do i = acoustic_surface(2,ispec_acoustic_surface), acoustic_surface(3,ispec_acoustic_surface)
@@ -46,11 +46,11 @@
potential_acoustic(iglob) = ZERO
potential_dot_acoustic(iglob) = ZERO
potential_dot_dot_acoustic(iglob) = ZERO
-
+
end do
end do
-
+
end do
-
+
end subroutine enforce_acoustic_free_surface
Modified: seismo/2D/SPECFEM2D/trunk/gll_library.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/gll_library.f90 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/gll_library.f90 2007-12-07 23:57:01 UTC (rev 8566)
@@ -487,7 +487,7 @@
integer np
double precision alpha,beta
- double precision z(np)
+ double precision z(np)
real(kind=CUSTOM_REAL) :: w(np)
integer n,nm1,i
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-07-27 16:08:31 UTC (rev 8565)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-07 23:57:01 UTC (rev 8566)
@@ -191,8 +191,11 @@
integer nspec_allocate
double precision :: deltatsquare,deltatcube,deltatfourth,twelvedeltat,fourdeltatsquare
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1,e11,e13
+ double precision, dimension(N_SLS) :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
+ double precision :: Mu_nu1,Mu_nu2
+
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
- e1_mech1,e11_mech1,e13_mech1,e1_mech2,e11_mech2,e13_mech2, &
dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1
! for fluid/solid coupling and edge detection
@@ -533,12 +536,10 @@
endif
! allocate memory variables for attenuation
- allocate(e1_mech1(NGLLX,NGLLZ,nspec_allocate))
- allocate(e11_mech1(NGLLX,NGLLZ,nspec_allocate))
- allocate(e13_mech1(NGLLX,NGLLZ,nspec_allocate))
- allocate(e1_mech2(NGLLX,NGLLZ,nspec_allocate))
- allocate(e11_mech2(NGLLX,NGLLZ,nspec_allocate))
- allocate(e13_mech2(NGLLX,NGLLZ,nspec_allocate))
+ allocate(e1(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(e11(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(e13(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+
allocate(dux_dxl_n(NGLLX,NGLLZ,nspec_allocate))
allocate(duz_dzl_n(NGLLX,NGLLZ,nspec_allocate))
allocate(duz_dxl_n(NGLLX,NGLLZ,nspec_allocate))
@@ -548,6 +549,8 @@
allocate(duz_dxl_np1(NGLLX,NGLLZ,nspec_allocate))
allocate(dux_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
+! define the attenuation constants
+ call attenuation_model(inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2)
!
!---- read interfaces data
@@ -1773,10 +1776,10 @@
initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,angleforce,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,elastic,codeabs, &
accel_elastic,veloc_elastic,displ_elastic,density,elastcoef,xix,xiz,gammax,gammaz, &
- jacobian,vpext,vsext,rhoext,source_time_function,sourcearray,e1_mech1,e11_mech1, &
- e13_mech1,e1_mech2,e11_mech2,e13_mech2,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
+ jacobian,vpext,vsext,rhoext,source_time_function,sourcearray, &
+ e1,e11,e13,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll)
+ hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2)
! *********************************************************
! ************* add coupling with the acoustic side
@@ -1889,8 +1892,8 @@
xix,xiz,gammax,gammaz,jacobian,ibool,elastic,hprime_xx,hprime_zz, &
nspec,npoin,assign_external_model,it,deltat,t0,kmato,elastcoef,density, &
vpext,vsext,rhoext,wxgll,wzgll,numat, &
- pressure_element,vector_field_element,e1_mech1,e11_mech1,e1_mech2,e11_mech2, &
- potential_dot_acoustic,potential_dot_dot_acoustic,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON)
+ pressure_element,vector_field_element,e1,e11, &
+ potential_dot_acoustic,potential_dot_dot_acoustic,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
!---- display time step and max of norm of displacement
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5) then
@@ -1931,8 +1934,8 @@
call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,elastic, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1_mech1,e11_mech1, &
- e1_mech2,e11_mech2,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON)
+ numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
+ TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
else if(.not. elastic(ispec)) then
@@ -2121,8 +2124,8 @@
call compute_pressure_whole_medium(potential_dot_dot_acoustic,displ_elastic,elastic,vector_field_display, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,elastcoef,vpext,vsext,rhoext,e1_mech1,e11_mech1, &
- e1_mech2,e11_mech2,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON)
+ numat,kmato,density,elastcoef,vpext,vsext,rhoext,e1,e11, &
+ TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
else
call exit_MPI('wrong type for snapshots')
More information about the cig-commits
mailing list