[cig-commits] r19804 - in seismo/2D/SPECFEM2D/trunk: DATA EXAMPLES EXAMPLES/Abel_Brest EXAMPLES/Gmsh_example_MPI EXAMPLES/Gmsh_example_serial EXAMPLES/INDUSTRIAL_FORMAT EXAMPLES/M2_UPPA EXAMPLES/Rayleigh_wave_no_crack EXAMPLES/Rayleigh_wave_with_crack EXAMPLES/Tape2007 EXAMPLES/Tape2007_kernel EXAMPLES/Tromp2005 EXAMPLES/Tromp2005_kernel EXAMPLES/acoustic_poroelastic EXAMPLES/attenuation EXAMPLES/canyon EXAMPLES/fluid_solid/fluid_solid_external_mesh EXAMPLES/fluid_solid/from_2000_Geophysics_paper_flat_ocean_bottom EXAMPLES/fluid_solid/from_2000_Geophysics_paper_sinusoidal_ocean_bottom EXAMPLES/init_plane EXAMPLES/noise_layered/model_0 EXAMPLES/noise_layered/model_1 EXAMPLES/noise_layered/model_2 EXAMPLES/noise_uniform EXAMPLES/semi_infinite_homo src/meshfem2D src/specfem2D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Mar 18 19:12:10 PDT 2012
Author: dkomati1
Date: 2012-03-18 19:12:09 -0700 (Sun, 18 Mar 2012)
New Revision: 19804
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Abel_Brest/Par_file_Abel_Balanche_bathy_source_solid
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_MPI/Par_file_Gmsh_SqrCircles.in
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/Par_file_Gmsh_SqrCircles.in
seismo/2D/SPECFEM2D/trunk/EXAMPLES/INDUSTRIAL_FORMAT/Par_file
seismo/2D/SPECFEM2D/trunk/EXAMPLES/M2_UPPA/Par_file_M2_UPPA
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_no_crack/Par_file_Rayleigh_2D
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/Par_file_Rayleigh_2D
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_132rec_checker
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_onerec
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/Par_file_Tape2007_onerec
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005_s100
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005_kernel/Par_file_Tromp2005
seismo/2D/SPECFEM2D/trunk/EXAMPLES/acoustic_poroelastic/Par_file_acoustic_poroelastic
seismo/2D/SPECFEM2D/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D
seismo/2D/SPECFEM2D/trunk/EXAMPLES/canyon/Par_file_canyon
seismo/2D/SPECFEM2D/trunk/EXAMPLES/change_all_2D_Par_files.csh
seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/fluid_solid_external_mesh/Par_file_fluid_solid
seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/from_2000_Geophysics_paper_flat_ocean_bottom/Par_file_fluid_solid
seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/from_2000_Geophysics_paper_sinusoidal_ocean_bottom/Par_file_fluid_solid
seismo/2D/SPECFEM2D/trunk/EXAMPLES/init_plane/Par_file_Slave
seismo/2D/SPECFEM2D/trunk/EXAMPLES/init_plane/Par_file_Slave_for
seismo/2D/SPECFEM2D/trunk/EXAMPLES/init_plane/Par_file_Slave_kernel
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_0/Par_file_fair
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_0/Par_file_good
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_1/Par_file_best
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_1/Par_file_fair
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_1/Par_file_good
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_2/Par_file_fair
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_2/Par_file_good
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3
seismo/2D/SPECFEM2D/trunk/EXAMPLES/semi_infinite_homo/Par_file_elastic_2D
seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_parameter_file.F90
seismo/2D/SPECFEM2D/trunk/src/meshfem2D/save_databases.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/assemble_MPI.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
added the missing C Delta_t / 2 term in the mass matrix for Clayton-Engquist absorbing elements (mentioned by Jean-Paul Ampuero in 2003, see the book of Hughes 1987 chapter 9).
Zhinan Xie also added a spring to further stabilize the Clayton-Engquist absorbing conditions.
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2012-03-19 02:12:09 UTC (rev 19804)
@@ -93,6 +93,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Abel_Brest/Par_file_Abel_Balanche_bathy_source_solid
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Abel_Brest/Par_file_Abel_Balanche_bathy_source_solid 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Abel_Brest/Par_file_Abel_Balanche_bathy_source_solid 2012-03-19 02:12:09 UTC (rev 19804)
@@ -99,6 +99,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_MPI/Par_file_Gmsh_SqrCircles.in
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_MPI/Par_file_Gmsh_SqrCircles.in 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_MPI/Par_file_Gmsh_SqrCircles.in 2012-03-19 02:12:09 UTC (rev 19804)
@@ -95,6 +95,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/Par_file_Gmsh_SqrCircles.in
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/Par_file_Gmsh_SqrCircles.in 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/Par_file_Gmsh_SqrCircles.in 2012-03-19 02:12:09 UTC (rev 19804)
@@ -95,6 +95,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/INDUSTRIAL_FORMAT/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/INDUSTRIAL_FORMAT/Par_file 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/INDUSTRIAL_FORMAT/Par_file 2012-03-19 02:12:09 UTC (rev 19804)
@@ -93,6 +93,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/M2_UPPA/Par_file_M2_UPPA
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/M2_UPPA/Par_file_M2_UPPA 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/M2_UPPA/Par_file_M2_UPPA 2012-03-19 02:12:09 UTC (rev 19804)
@@ -93,6 +93,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_no_crack/Par_file_Rayleigh_2D
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_no_crack/Par_file_Rayleigh_2D 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_no_crack/Par_file_Rayleigh_2D 2012-03-19 02:12:09 UTC (rev 19804)
@@ -90,6 +90,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/Par_file_Rayleigh_2D
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/Par_file_Rayleigh_2D 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/Par_file_Rayleigh_2D 2012-03-19 02:12:09 UTC (rev 19804)
@@ -94,6 +94,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_132rec_checker
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_132rec_checker 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_132rec_checker 2012-03-19 02:12:09 UTC (rev 19804)
@@ -876,6 +876,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_onerec
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_onerec 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_onerec 2012-03-19 02:12:09 UTC (rev 19804)
@@ -91,6 +91,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/Par_file_Tape2007_onerec
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/Par_file_Tape2007_onerec 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/Par_file_Tape2007_onerec 2012-03-19 02:12:09 UTC (rev 19804)
@@ -91,6 +91,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005 2012-03-19 02:12:09 UTC (rev 19804)
@@ -90,6 +90,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005_s100
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005_s100 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005_s100 2012-03-19 02:12:09 UTC (rev 19804)
@@ -90,6 +90,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005_kernel/Par_file_Tromp2005
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005_kernel/Par_file_Tromp2005 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005_kernel/Par_file_Tromp2005 2012-03-19 02:12:09 UTC (rev 19804)
@@ -90,6 +90,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/acoustic_poroelastic/Par_file_acoustic_poroelastic
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/acoustic_poroelastic/Par_file_acoustic_poroelastic 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/acoustic_poroelastic/Par_file_acoustic_poroelastic 2012-03-19 02:12:09 UTC (rev 19804)
@@ -99,6 +99,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D 2012-03-19 02:12:09 UTC (rev 19804)
@@ -90,6 +90,7 @@
# absorbing boundary active or not
absorbing_conditions = .false.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/canyon/Par_file_canyon
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/canyon/Par_file_canyon 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/canyon/Par_file_canyon 2012-03-19 02:12:09 UTC (rev 19804)
@@ -90,6 +90,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/change_all_2D_Par_files.csh
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/change_all_2D_Par_files.csh 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/change_all_2D_Par_files.csh 2012-03-19 02:12:09 UTC (rev 19804)
@@ -11,11 +11,12 @@
set newfile = "../"$file
# DK DK leave the white spaces, in order to have the right alignment with other variables
+###sed -e '1,$s/PERFORM_CUTHILL_MCKEE = .true. # perform inverse Cuthill-McKee (1969) optimization\/permutation for mesh numbering/PERFORM_CUTHILL_MCKEE = .false. # perform inverse Cuthill-McKee (1969) optimization\/permutation for mesh numbering (can be very costly and not very useful)/g' < $newfile > __________temp_27_zzzyyyyxxxx__________
sed -e '1,$s/PERFORM_CUTHILL_MCKEE = .true. # perform inverse Cuthill-McKee (1969) optimization\/permutation for mesh numbering/PERFORM_CUTHILL_MCKEE = .false. # perform inverse Cuthill-McKee (1969) optimization\/permutation for mesh numbering (can be very costly and not very useful)/g' < $newfile > __________temp_27_zzzyyyyxxxx__________
clean_a_Par_file.py __________temp_27_zzzyyyyxxxx__________ $file
-mv __________temp_27_zzzyyyyxxxx__________ $newfile
+#mv __________temp_27_zzzyyyyxxxx__________ $newfile
end
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/fluid_solid_external_mesh/Par_file_fluid_solid
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/fluid_solid_external_mesh/Par_file_fluid_solid 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/fluid_solid_external_mesh/Par_file_fluid_solid 2012-03-19 02:12:09 UTC (rev 19804)
@@ -95,6 +95,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/from_2000_Geophysics_paper_flat_ocean_bottom/Par_file_fluid_solid
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/from_2000_Geophysics_paper_flat_ocean_bottom/Par_file_fluid_solid 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/from_2000_Geophysics_paper_flat_ocean_bottom/Par_file_fluid_solid 2012-03-19 02:12:09 UTC (rev 19804)
@@ -92,6 +92,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/from_2000_Geophysics_paper_sinusoidal_ocean_bottom/Par_file_fluid_solid
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/from_2000_Geophysics_paper_sinusoidal_ocean_bottom/Par_file_fluid_solid 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/from_2000_Geophysics_paper_sinusoidal_ocean_bottom/Par_file_fluid_solid 2012-03-19 02:12:09 UTC (rev 19804)
@@ -92,6 +92,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/init_plane/Par_file_Slave
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/init_plane/Par_file_Slave 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/init_plane/Par_file_Slave 2012-03-19 02:12:09 UTC (rev 19804)
@@ -91,6 +91,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/init_plane/Par_file_Slave_for
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/init_plane/Par_file_Slave_for 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/init_plane/Par_file_Slave_for 2012-03-19 02:12:09 UTC (rev 19804)
@@ -91,6 +91,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/init_plane/Par_file_Slave_kernel
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/init_plane/Par_file_Slave_kernel 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/init_plane/Par_file_Slave_kernel 2012-03-19 02:12:09 UTC (rev 19804)
@@ -91,6 +91,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_0/Par_file_fair
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_0/Par_file_fair 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_0/Par_file_fair 2012-03-19 02:12:09 UTC (rev 19804)
@@ -90,6 +90,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_0/Par_file_good
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_0/Par_file_good 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_0/Par_file_good 2012-03-19 02:12:09 UTC (rev 19804)
@@ -90,6 +90,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_1/Par_file_best
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_1/Par_file_best 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_1/Par_file_best 2012-03-19 02:12:09 UTC (rev 19804)
@@ -91,6 +91,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_1/Par_file_fair
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_1/Par_file_fair 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_1/Par_file_fair 2012-03-19 02:12:09 UTC (rev 19804)
@@ -91,6 +91,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_1/Par_file_good
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_1/Par_file_good 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_1/Par_file_good 2012-03-19 02:12:09 UTC (rev 19804)
@@ -91,6 +91,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_2/Par_file_fair
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_2/Par_file_fair 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_2/Par_file_fair 2012-03-19 02:12:09 UTC (rev 19804)
@@ -92,6 +92,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_2/Par_file_good
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_2/Par_file_good 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_layered/model_2/Par_file_good 2012-03-19 02:12:09 UTC (rev 19804)
@@ -92,6 +92,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1 2012-03-19 02:12:09 UTC (rev 19804)
@@ -90,6 +90,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2 2012-03-19 02:12:09 UTC (rev 19804)
@@ -90,6 +90,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3 2012-03-19 02:12:09 UTC (rev 19804)
@@ -90,6 +90,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/semi_infinite_homo/Par_file_elastic_2D
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/semi_infinite_homo/Par_file_elastic_2D 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/semi_infinite_homo/Par_file_elastic_2D 2012-03-19 02:12:09 UTC (rev 19804)
@@ -90,6 +90,7 @@
# absorbing boundary active or not
absorbing_conditions = .true.
+ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
Modified: seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_parameter_file.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_parameter_file.F90 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_parameter_file.F90 2012-03-19 02:12:09 UTC (rev 19804)
@@ -74,6 +74,7 @@
logical :: p_sv
logical :: any_abs,absbottom,absright,abstop,absleft
+ logical :: ADD_SPRING_TO_STACEY
integer :: nt
double precision :: deltat
@@ -239,8 +240,8 @@
call read_value_double_precision_p(USER_T0, 'solver.USER_T0')
if(err_occurred() /= 0) stop 'error reading parameter 17b in Par_file'
- call read_value_integer_p(time_stepping_scheme, 'solver.time_stepping_scheme') !xiezhinan
- if(err_occurred() /= 0) stop 'error reading parameter 17c in Par_file' !xiezhinan
+ call read_value_integer_p(time_stepping_scheme, 'solver.time_stepping_scheme')
+ if(err_occurred() /= 0) stop 'error reading parameter 17c in Par_file'
! read source infos
call read_value_integer_p(NSOURCES, 'solver.NSOURCES')
@@ -419,6 +420,9 @@
call read_value_logical_p(any_abs, 'solver.absorbing_conditions')
if(err_occurred() /= 0) stop 'error reading parameter 51a in Par_file'
+ call read_value_logical_p(ADD_SPRING_TO_STACEY, 'solver.ADD_SPRING_TO_STACEY')
+ if(err_occurred() /= 0) stop 'error reading parameter 51a in Par_file'
+
call read_value_logical_p(ADD_PERIODIC_CONDITIONS, 'solver.ADD_PERIODIC_CONDITIONS')
if(err_occurred() /= 0) stop 'error reading parameter 51b in Par_file'
Modified: seismo/2D/SPECFEM2D/trunk/src/meshfem2D/save_databases.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/meshfem2D/save_databases.f90 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/src/meshfem2D/save_databases.f90 2012-03-19 02:12:09 UTC (rev 19804)
@@ -166,6 +166,9 @@
write(15,*) 'time_stepping_scheme'
write(15,*) time_stepping_scheme
+ write(15,*) 'ADD_SPRING_TO_STACEY'
+ write(15,*) ADD_SPRING_TO_STACEY
+
write(15,*) 'ADD_PERIODIC_CONDITIONS'
write(15,*) ADD_PERIODIC_CONDITIONS
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/assemble_MPI.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/assemble_MPI.F90 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/assemble_MPI.F90 2012-03-19 02:12:09 UTC (rev 19804)
@@ -55,7 +55,7 @@
! Assembling the mass matrix.
!-----------------------------------------------
subroutine assemble_MPI_scalar(array_val1,npoin_val1, &
- array_val2,npoin_val2, &
+ array_val2,array_val5,npoin_val2, &
array_val3,array_val4,npoin_val3, &
ninterface, max_interface_size, max_ibool_interfaces_size_ac, &
max_ibool_interfaces_size_el, &
@@ -85,7 +85,7 @@
real(kind=CUSTOM_REAL), dimension(npoin_val1), intent(inout) :: array_val1
! elastic
integer :: npoin_val2
- real(kind=CUSTOM_REAL), dimension(npoin_val2), intent(inout) :: array_val2
+ real(kind=CUSTOM_REAL), dimension(npoin_val2), intent(inout) :: array_val2,array_val5
! poroelastic
integer :: npoin_val3
real(kind=CUSTOM_REAL), dimension(npoin_val3), intent(inout) :: array_val3,array_val4
@@ -118,6 +118,12 @@
array_val2(ibool_interfaces_elastic(i,num_interface))
end do
+ do i = 1, nibool_interfaces_elastic(num_interface)
+ ipoin = ipoin + 1
+ buffer_send_faces_scalar(ipoin,num_interface) = &
+ array_val5(ibool_interfaces_elastic(i,num_interface))
+ end do
+
do i = 1, nibool_interfaces_poroelastic(num_interface)
ipoin = ipoin + 1
buffer_send_faces_scalar(ipoin,num_interface) = &
@@ -164,6 +170,13 @@
+ buffer_recv_faces_scalar(ipoin,num_interface)
end do
+ do i = 1, nibool_interfaces_elastic(num_interface)
+ ipoin = ipoin + 1
+ array_val5(ibool_interfaces_elastic(i,num_interface)) = &
+ array_val5(ibool_interfaces_elastic(i,num_interface)) &
+ + buffer_recv_faces_scalar(ipoin,num_interface)
+ end do
+
do i = 1, nibool_interfaces_poroelastic(num_interface)
ipoin = ipoin + 1
array_val3(ibool_interfaces_poroelastic(i,num_interface)) = &
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90 2012-03-19 02:12:09 UTC (rev 19804)
@@ -62,7 +62,7 @@
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,mu_k,kappa_k,&
e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_RK, e11_force_RK, e13_force_RK, &
- stage_time_scheme,i_stage)
+ stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring)
! compute forces for the elastic elements
@@ -86,6 +86,8 @@
integer :: stage_time_scheme,i_stage
logical :: anyabs,assign_external_model,initialfield,ATTENUATION_VISCOELASTIC_SOLID,add_Bielak_conditions
+ logical :: ADD_SPRING_TO_STACEY
+ real(kind=CUSTOM_REAL) :: x_center_spring,z_center_spring
logical :: SAVE_FORWARD
@@ -163,6 +165,7 @@
real(kind=CUSTOM_REAL) :: sigma_xx,sigma_xy,sigma_xz,sigma_zy,sigma_zz,sigma_zx
real(kind=CUSTOM_REAL) :: b_sigma_xx,b_sigma_xy,b_sigma_xz,b_sigma_zy,b_sigma_zz,b_sigma_zx
real(kind=CUSTOM_REAL) :: nx,nz,vx,vy,vz,vn,rho_vp,rho_vs,tx,ty,tz,weight,xxi,zxi,xgamma,zgamma,jacobian1D
+ real(kind=CUSTOM_REAL) :: displx,disply,displz,displn,spring_position,displtx,displty,displtz
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: tempx1,tempx2,tempy1,tempy2,tempz1,tempz2
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: b_tempx1,b_tempx2,b_tempy1,b_tempy2,b_tempz1,b_tempz2
@@ -577,9 +580,34 @@
ty = rho_vs*vy
tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0)*weight
+ displtx=0.0d0
+ displtz=0.0d0
+
+ if(ADD_SPRING_TO_STACEY)then
+
+ displx = displ_elastic(1,iglob)
+ disply = displ_elastic(2,iglob)
+ displz = displ_elastic(3,iglob)
+
+ spring_position=sqrt((coord(1,iglob)-x_center_spring)**2 +&
+ (coord(2,iglob)-z_center_spring)**2)
+
+ displn = nx*displx+nz*displz
+
+ displtx = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)/&
+ (2.0*spring_position)*displn*nx &
+ +mul_unrelaxed_elastic/(2.0*spring_position)*(displx-displn*nx)
+ displty = mul_unrelaxed_elastic*disply
+ displtz = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)/&
+ (2.0*spring_position)*displn*nz &
+ +mul_unrelaxed_elastic/(2.0*spring_position)*(displz-displn*nz)
+
+ endif
+
+
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0+displtx)*weight
accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
- accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0)*weight
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0+displtz)*weight
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
if(p_sv)then !P-SV waves
@@ -668,9 +696,34 @@
ty = rho_vs*vy
tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx - traction_x_t0)*weight
+ displtx=0.0d0
+ displtz=0.0d0
+
+ if(ADD_SPRING_TO_STACEY)then
+
+ displx = displ_elastic(1,iglob)
+ disply = displ_elastic(2,iglob)
+ displz = displ_elastic(3,iglob)
+
+ spring_position=sqrt((coord(1,iglob)-x_center_spring)**2 +&
+ (coord(2,iglob)-z_center_spring)**2)
+
+ displn = nx*displx+nz*displz
+
+ displtx = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)/&
+ (2.0*spring_position)*displn*nx &
+ +mul_unrelaxed_elastic/(2.0*spring_position)*(displx-displn*nx)
+ displty = mul_unrelaxed_elastic*disply
+ displtz = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)/&
+ (2.0*spring_position)*displn*nz &
+ +mul_unrelaxed_elastic/(2.0*spring_position)*(displz-displn*nz)
+
+ endif
+
+
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0+displtx)*weight
accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
- accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz - traction_z_t0)*weight
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0+displtz)*weight
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
if(p_sv)then !P-SV waves
@@ -774,9 +827,41 @@
tz = 0
endif
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0)*weight
+ displtx=0.0d0
+ displtz=0.0d0
+
+ if(ADD_SPRING_TO_STACEY)then
+
+ displx = displ_elastic(1,iglob)
+ disply = displ_elastic(2,iglob)
+ displz = displ_elastic(3,iglob)
+
+ spring_position=sqrt((coord(1,iglob)-x_center_spring)**2 +&
+ (coord(2,iglob)-z_center_spring)**2)
+
+ displn = nx*displx+nz*displz
+
+ displtx = (lambdal_unrelaxed_elastic+2.0*mul_unrelaxed_elastic)/&
+ (2.0*spring_position)*displn*nx &
+ +mul_unrelaxed_elastic/(2.0*spring_position)*(displx-displn*nx)
+ displty = mul_unrelaxed_elastic*disply
+ displtz = (lambdal_unrelaxed_elastic+2.0*mul_unrelaxed_elastic)/&
+ (2.0*spring_position)*displn*nz &
+ +mul_unrelaxed_elastic/(2.0*spring_position)*(displz-displn*nz)
+
+ if((codeabs(ILEFT,ispecabs) .and. i == 1) .or. (codeabs(IRIGHT,ispecabs) .and. i == NGLLX)) then
+ displtx = 0
+ displty = 0
+ displtz = 0
+ endif
+
+
+ endif
+
+
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0+displtx)*weight
accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
- accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0)*weight
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0+displtz)*weight
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
if(p_sv)then !P-SV waves
@@ -872,9 +957,41 @@
tz = 0
endif
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx - traction_x_t0)*weight
+ displtx=0.0d0
+ displtz=0.0d0
+
+ if(ADD_SPRING_TO_STACEY)then
+
+ displx = displ_elastic(1,iglob)
+ disply = displ_elastic(2,iglob)
+ displz = displ_elastic(3,iglob)
+
+ spring_position=sqrt((coord(1,iglob)-x_center_spring)**2 +&
+ (coord(2,iglob)-z_center_spring)**2)
+
+ displn = nx*displx+nz*displz
+
+ displtx = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)/&
+ (2.0*spring_position)*displn*nx &
+ +mul_unrelaxed_elastic/(2.0*spring_position)*(displx-displn*nx)
+ displty = mul_unrelaxed_elastic*disply
+ displtz = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)/&
+ (2.0*spring_position)*displn*nz &
+ +mul_unrelaxed_elastic/(2.0*spring_position)*(displz-displn*nz)
+
+ if((codeabs(ILEFT,ispecabs) .and. i == 1) .or. (codeabs(IRIGHT,ispecabs) .and. i == NGLLX)) then
+ displtx = 0
+ displty = 0
+ displtz = 0
+ endif
+
+
+ endif
+
+
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0+displtx)*weight
accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
- accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz - traction_z_t0)*weight
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0+displtz)*weight
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
if(p_sv)then !P-SV waves
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.f90 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.f90 2012-03-19 02:12:09 UTC (rev 19804)
@@ -44,7 +44,7 @@
!========================================================================
subroutine invert_mass_matrix_init(any_elastic,any_acoustic,any_poroelastic, &
- rmass_inverse_elastic,nglob_elastic, &
+ rmass_inverse_elastic_one,nglob_elastic, &
rmass_inverse_acoustic,nglob_acoustic, &
rmass_s_inverse_poroelastic, &
rmass_w_inverse_poroelastic,nglob_poroelastic, &
@@ -52,19 +52,40 @@
elastic,poroelastic, &
assign_external_model,numat, &
density,poroelastcoef,porosity,tortuosity, &
- vpext,rhoext)
+ vpext,rhoext,&
+ anyabs,numabs,deltat,codeabs,rmass_inverse_elastic_three,&
+ nelemabs,vsext,xix,xiz,gammaz,gammax)
! builds the global mass matrix
implicit none
include 'constants.h'
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ logical :: anyabs
+ integer :: nelemabs,ibegin,iend,ispecabs
+! integer :: ispec,i,j,k,iglob,ispecabs,ibegin,iend,irec,irec_local
+ integer, dimension(nelemabs) :: numabs
+ double precision :: deltat
+ logical, dimension(4,nelemabs) :: codeabs
+
+!!local parameter
+ ! material properties of the elastic medium
+ real(kind=CUSTOM_REAL) :: mul_unrelaxed_elastic,lambdal_unrelaxed_elastic,cpl,csl
+ integer count_left,count_right,count_bottom
+ real(kind=CUSTOM_REAL) :: nx,nz,vx,vy,vz,vn,rho_vp,rho_vs,tx,ty,tz,&
+ weight,xxi,zxi,xgamma,zgamma,jacobian1D
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
logical any_elastic,any_acoustic,any_poroelastic
! inverse mass matrices
integer :: nglob_elastic
- real(kind=CUSTOM_REAL), dimension(nglob_elastic) :: rmass_inverse_elastic
+ real(kind=CUSTOM_REAL), dimension(nglob_elastic) :: rmass_inverse_elastic_one,& !zhinan
+ rmass_inverse_elastic_three
+
integer :: nglob_acoustic
real(kind=CUSTOM_REAL), dimension(nglob_acoustic) :: rmass_inverse_acoustic
@@ -86,15 +107,18 @@
double precision, dimension(2,numat) :: density
double precision, dimension(4,3,numat) :: poroelastcoef
double precision, dimension(numat) :: porosity,tortuosity
- double precision, dimension(NGLLX,NGLLX,nspec) :: vpext,rhoext
+ double precision, dimension(NGLLX,NGLLX,nspec) :: vpext,rhoext,vsext !zhinan
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz !zhinan
+
! local parameters
integer :: ispec,i,j,iglob
double precision :: rhol,kappal,mul_relaxed,lambdal_relaxed
double precision :: rhol_s,rhol_f,rhol_bar,phil,tortl
! initializes mass matrix
- if(any_elastic) rmass_inverse_elastic(:) = 0._CUSTOM_REAL
+ if(any_elastic) rmass_inverse_elastic_one(:) = 0._CUSTOM_REAL
+ if(any_elastic) rmass_inverse_elastic_three(:) = 0._CUSTOM_REAL
if(any_poroelastic) rmass_s_inverse_poroelastic(:) = 0._CUSTOM_REAL
if(any_poroelastic) rmass_w_inverse_poroelastic(:) = 0._CUSTOM_REAL
if(any_acoustic) rmass_inverse_acoustic(:) = 0._CUSTOM_REAL
@@ -136,10 +160,11 @@
elseif( elastic(ispec) ) then
! for elastic medium
-
- rmass_inverse_elastic(iglob) = rmass_inverse_elastic(iglob) &
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec)
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+
else
! for acoustic medium
@@ -153,13 +178,270 @@
enddo
enddo ! do ispec = 1,nspec
+ !
+ !--- DK and Zhinan Xie: add C Delta_t / 2 contribution to the mass matrix
+ !--- DK and Zhinan Xie: in the case of Clayton-Engquist absorbing boundaries;
+ !--- DK and Zhinan Xie: see for instance the book of Hughes (1987) chapter 9.
+ !--- DK and Zhinan Xie: IMPORTANT: note that this implies that we must have two different mass matrices,
+ !--- DK and Zhinan Xie: one per component of the wave field i.e. one per spatial dimension.
+ !--- DK and Zhinan Xie: This was also suggested by Jean-Paul Ampuero in 2003.
+ !
+ if(anyabs) then
+ count_left=1
+ count_right=1
+ count_bottom=1
+ do ispecabs = 1,nelemabs
+ ispec = numabs(ispecabs)
+ ! get elastic parameters of current spectral elemegammaznt
+ lambdal_unrelaxed_elastic = poroelastcoef(1,1,kmato(ispec))
+ mul_unrelaxed_elastic = poroelastcoef(2,1,kmato(ispec))
+ rhol = density(1,kmato(ispec))
+ kappal = lambdal_unrelaxed_elastic + TWO*mul_unrelaxed_elastic/3._CUSTOM_REAL
+ cpl = sqrt((kappal + 4._CUSTOM_REAL*mul_unrelaxed_elastic/3._CUSTOM_REAL)/rhol)
+ csl = sqrt(mul_unrelaxed_elastic/rhol)
+
+ !--- left absorbing boundary
+ if(codeabs(ILEFT,ispecabs)) then
+
+ i = 1
+
+ do j = 1,NGLLZ
+
+ iglob = ibool(i,j,ispec)
+
+ ! external velocity model
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ csl = vsext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
+
+ rho_vp = rhol*cpl
+ rho_vs = rhol*csl
+
+ xgamma = - xiz(i,j,ispec) * jacobian(i,j,ispec)
+ zgamma = + xix(i,j,ispec) * jacobian(i,j,ispec)
+ jacobian1D = sqrt(xgamma**2 + zgamma**2)
+ nx = - zgamma / jacobian1D
+ nz = + xgamma / jacobian1D
+
+ weight = jacobian1D * wzgll(j)
+
+ ! Clayton-Engquist condition if elastic
+ if(elastic(ispec)) then
+
+ vx = 1.0d0*deltat/2.0d0
+ vy = 1.0d0*deltat/2.0d0
+ vz = 1.0d0*deltat/2.0d0
+
+ vn = nx*vx+nz*vz
+
+ tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ ty = rho_vs*vy
+ tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
+
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + (tx)*weight
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob) &
+ + (tz)*weight
+
+ endif
+ enddo
+
+ endif ! end of left absorbing boundary
+
+ !--- right absorbing boundary
+ if(codeabs(IRIGHT,ispecabs)) then
+
+ i = NGLLX
+
+ do j = 1,NGLLZ
+
+ iglob = ibool(i,j,ispec)
+
+ ! for analytical initial plane wave for Bielak's conditions
+ ! left or right edge, horizontal normal vector
+
+ ! external velocity model
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ csl = vsext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
+
+ rho_vp = rhol*cpl
+ rho_vs = rhol*csl
+
+ xgamma = - xiz(i,j,ispec) * jacobian(i,j,ispec)
+ zgamma = + xix(i,j,ispec) * jacobian(i,j,ispec)
+ jacobian1D = sqrt(xgamma**2 + zgamma**2)
+ nx = + zgamma / jacobian1D
+ nz = - xgamma / jacobian1D
+
+ weight = jacobian1D * wzgll(j)
+
+ ! Clayton-Engquist condition if elastic
+ if(elastic(ispec)) then
+
+ vx = 1.0d0*deltat/2.0d0
+ vy = 1.0d0*deltat/2.0d0
+ vz = 1.0d0*deltat/2.0d0
+
+ vn = nx*vx+nz*vz
+
+ tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ ty = rho_vs*vy
+ tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
+
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + (tx)*weight
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob) &
+ + (tz)*weight
+
+ endif
+
+ enddo
+
+ endif ! end of right absorbing boundary
+
+ !--- bottom absorbing boundary
+ if(codeabs(IBOTTOM,ispecabs)) then
+
+ j = 1
+ ibegin = 1
+ iend = NGLLX
+
+ do i = ibegin,iend
+
+ iglob = ibool(i,j,ispec)
+ ! external velocity model
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ csl = vsext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
+
+ rho_vp = rhol*cpl
+ rho_vs = rhol*csl
+
+ xxi = + gammaz(i,j,ispec) * jacobian(i,j,ispec)
+ zxi = - gammax(i,j,ispec) * jacobian(i,j,ispec)
+ jacobian1D = sqrt(xxi**2 + zxi**2)
+ nx = + zxi / jacobian1D
+ nz = - xxi / jacobian1D
+
+ weight = jacobian1D * wxgll(i)
+
+ ! Clayton-Engquist condition if elastic
+ if(elastic(ispec)) then
+
+ vx = 1.0d0*deltat/2.0d0
+ vy = 1.0d0*deltat/2.0d0
+ vz = 1.0d0*deltat/2.0d0
+
+ vn = nx*vx+nz*vz
+
+ tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ ty = rho_vs*vy
+ tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
+
+! exclude corners to make sure there is no contradiction on the normal
+! for Stacey absorbing conditions but not for incident plane waves;
+! thus subtract nothing i.e. zero in that case
+ if((codeabs(ILEFT,ispecabs) .and. i == 1) .or. (codeabs(IRIGHT,ispecabs) .and. i == NGLLX)) then
+ tx = 0
+ ty = 0
+ tz = 0
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob)
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_three(iglob)
+ else
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + (tx)*weight
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob) &
+ + (tz)*weight
+ endif
+
+ endif
+
+ enddo
+
+ endif ! end of bottom absorbing boundary
+
+ !--- top absorbing boundary
+ if(codeabs(ITOP,ispecabs)) then
+
+ j = NGLLZ
+
+ ibegin = 1
+ iend = NGLLX
+
+ do i = ibegin,iend
+
+ iglob = ibool(i,j,ispec)
+
+ ! external velocity model
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ csl = vsext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
+
+ rho_vp = rhol*cpl
+ rho_vs = rhol*csl
+
+ xxi = + gammaz(i,j,ispec) * jacobian(i,j,ispec)
+ zxi = - gammax(i,j,ispec) * jacobian(i,j,ispec)
+ jacobian1D = sqrt(xxi**2 + zxi**2)
+ nx = - zxi / jacobian1D
+ nz = + xxi / jacobian1D
+
+ weight = jacobian1D * wxgll(i)
+
+ ! Clayton-Engquist condition if elastic
+ if(elastic(ispec)) then
+
+ vx = 1.0d0*deltat/2.0d0
+ vy = 1.0d0*deltat/2.0d0
+ vz = 1.0d0*deltat/2.0d0
+
+ vn = nx*vx+nz*vz
+
+ tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ ty = rho_vs*vy
+ tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
+
+! exclude corners to make sure there is no contradiction on the normal
+! for Stacey absorbing conditions but not for incident plane waves;
+! thus subtract nothing i.e. zero in that case
+ if((codeabs(ILEFT,ispecabs) .and. i == 1) .or. (codeabs(IRIGHT,ispecabs) .and. i == NGLLX)) then
+ tx = 0
+ ty = 0
+ tz = 0
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob)
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_three(iglob)
+ else
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + (tx)*weight
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob) &
+ + (tz)*weight
+ endif
+ endif
+
+
+ enddo
+ endif ! end of top absorbing boundary
+ enddo
+ endif ! end of absorbing boundaries
+
+
end subroutine invert_mass_matrix_init
!
!-------------------------------------------------------------------------------------------------
!
subroutine invert_mass_matrix(any_elastic,any_acoustic,any_poroelastic, &
- rmass_inverse_elastic,nglob_elastic, &
+ rmass_inverse_elastic_one,rmass_inverse_elastic_three,&
+ nglob_elastic, &
rmass_inverse_acoustic,nglob_acoustic, &
rmass_s_inverse_poroelastic, &
rmass_w_inverse_poroelastic,nglob_poroelastic)
@@ -173,7 +455,8 @@
! inverse mass matrices
integer :: nglob_elastic
- real(kind=CUSTOM_REAL), dimension(nglob_elastic) :: rmass_inverse_elastic
+ real(kind=CUSTOM_REAL), dimension(nglob_elastic) :: rmass_inverse_elastic_one,&
+ rmass_inverse_elastic_three
integer :: nglob_acoustic
real(kind=CUSTOM_REAL), dimension(nglob_acoustic) :: rmass_inverse_acoustic
@@ -185,7 +468,9 @@
! fill mass matrix with fictitious non-zero values to make sure it can be inverted globally
if(any_elastic) &
- where(rmass_inverse_elastic <= 0._CUSTOM_REAL) rmass_inverse_elastic = 1._CUSTOM_REAL
+ where(rmass_inverse_elastic_one <= 0._CUSTOM_REAL) rmass_inverse_elastic_one = 1._CUSTOM_REAL
+ if(any_elastic) &
+ where(rmass_inverse_elastic_three <= 0._CUSTOM_REAL) rmass_inverse_elastic_three = 1._CUSTOM_REAL
if(any_poroelastic) &
where(rmass_s_inverse_poroelastic <= 0._CUSTOM_REAL) rmass_s_inverse_poroelastic = 1._CUSTOM_REAL
if(any_poroelastic) &
@@ -195,7 +480,9 @@
! compute the inverse of the mass matrix
if(any_elastic) &
- rmass_inverse_elastic(:) = 1._CUSTOM_REAL / rmass_inverse_elastic(:)
+ rmass_inverse_elastic_one(:) = 1._CUSTOM_REAL / rmass_inverse_elastic_one(:)
+ if(any_elastic) &
+ rmass_inverse_elastic_three(:) = 1._CUSTOM_REAL / rmass_inverse_elastic_three(:)
if(any_poroelastic) &
rmass_s_inverse_poroelastic(:) = 1._CUSTOM_REAL / rmass_s_inverse_poroelastic(:)
if(any_poroelastic) &
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90 2012-03-19 02:12:09 UTC (rev 19804)
@@ -56,7 +56,7 @@
NSTEP,deltat,NTSTEP_BETWEEN_OUTPUT_SEISMO,NSOURCES, &
factor_subsample_image,USE_SNAPSHOT_NUMBER_IN_FILENAME,DRAW_WATER_CONSTANT_BLUE_IN_JPG,US_LETTER, &
POWER_DISPLAY_COLOR,PERFORM_CUTHILL_MCKEE,SU_FORMAT,USER_T0,time_stepping_scheme,&
- ADD_PERIODIC_CONDITIONS,PERIODIC_horiz_dist,PERIODIC_DETECT_TOL)
+ ADD_SPRING_TO_STACEY,ADD_PERIODIC_CONDITIONS,PERIODIC_horiz_dist,PERIODIC_DETECT_TOL)
! starts reading in parameters from input Database file
@@ -114,6 +114,9 @@
! 3 = classical 4th-order 4-stage Runge-Kutta
integer :: time_stepping_scheme
+!! DK DK for add spring to stacey absorbing boundary condition
+ logical :: ADD_SPRING_TO_STACEY
+
!! DK DK for horizontal periodic conditions: detect common points between left and right edges
logical :: ADD_PERIODIC_CONDITIONS
@@ -238,6 +241,9 @@
read(IIN,*) time_stepping_scheme
read(IIN,"(a80)") datlin
+ read(IIN,*) ADD_SPRING_TO_STACEY
+
+ read(IIN,"(a80)") datlin
read(IIN,*) ADD_PERIODIC_CONDITIONS
read(IIN,"(a80)") datlin
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-03-19 01:25:41 UTC (rev 19803)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-03-19 02:12:09 UTC (rev 19804)
@@ -482,7 +482,9 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: potential_acoustic_adj_coupling
! inverse mass matrices
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_elastic
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_elastic_one !zhinan
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_elastic_three !zhinan
+
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_acoustic
real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic
@@ -838,6 +840,10 @@
! to help locate elements with a negative Jacobian using OpenDX
logical :: found_a_negative_jacobian
+!! DK DK for add spring to stacey absorbing boundary condition
+ logical :: ADD_SPRING_TO_STACEY
+ double precision :: x_center_spring,z_center_spring,xmin_spring,xmax_spring,zmin_spring,zmax_spring
+
!! DK DK for horizontal periodic conditions: detect common points between left and right edges
logical :: ADD_PERIODIC_CONDITIONS
@@ -967,7 +973,7 @@
NSTEP,deltat,NTSTEP_BETWEEN_OUTPUT_SEISMO,NSOURCES, &
factor_subsample_image,USE_SNAPSHOT_NUMBER_IN_FILENAME,DRAW_WATER_CONSTANT_BLUE_IN_JPG,US_LETTER, &
POWER_DISPLAY_COLOR,PERFORM_CUTHILL_MCKEE,SU_FORMAT,USER_T0, time_stepping_scheme, &
- ADD_PERIODIC_CONDITIONS,PERIODIC_horiz_dist,PERIODIC_DETECT_TOL)
+ ADD_SPRING_TO_STACEY,ADD_PERIODIC_CONDITIONS,PERIODIC_horiz_dist,PERIODIC_DETECT_TOL)
if(nproc_read_from_database < 1) stop 'should have nproc_read_from_database >= 1'
if(SIMULATION_TYPE == 2 .and.(time_stepping_scheme == 2 .or. time_stepping_scheme == 3)) &
stop 'RK and LDDRK time scheme not supported for adjoint inversion'
@@ -1001,7 +1007,7 @@
NSTEP,deltat,NTSTEP_BETWEEN_OUTPUT_SEISMO,NSOURCES, &
factor_subsample_image,USE_SNAPSHOT_NUMBER_IN_FILENAME,DRAW_WATER_CONSTANT_BLUE_IN_JPG,US_LETTER, &
POWER_DISPLAY_COLOR,PERFORM_CUTHILL_MCKEE,SU_FORMAT,USER_T0, time_stepping_scheme, &
- ADD_PERIODIC_CONDITIONS,PERIODIC_horiz_dist,PERIODIC_DETECT_TOL)
+ ADD_SPRING_TO_STACEY,ADD_PERIODIC_CONDITIONS,PERIODIC_horiz_dist,PERIODIC_DETECT_TOL)
!
!--- source information
@@ -1612,6 +1618,13 @@
enddo
enddo
+ xmin_spring = minval(coord(1,:))
+ xmax_spring = maxval(coord(1,:))
+ zmin_spring = minval(coord(2,:))
+ zmax_spring = maxval(coord(2,:))
+
+ x_center_spring=(xmax_spring+xmin_spring)/2.d0
+ z_center_spring=(zmax_spring+zmin_spring)/2.d0
!
!---- generate the global numbering
!
@@ -1741,7 +1754,7 @@
!
!---- build the global mass matrix and invert it once and for all
!
- rmass_inverse_elastic(:) = 0._CUSTOM_REAL
+ rmass_inverse_elastic_one(:) = 0._CUSTOM_REAL
do ispec = 1,nspec
do j = 1,NGLLZ
do i = 1,NGLLX
@@ -1758,7 +1771,7 @@
kappal = lambdal_unrelaxed_elastic + 2.d0/3.d0*mul_unrelaxed_elastic
endif
- rmass_inverse_elastic(iglob) = rmass_inverse_elastic(iglob) &
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec)
enddo
@@ -1769,8 +1782,8 @@
! set entries that are equal to zero to something else, e.g. 1, to avoid division by zero
! these degrees of freedom correspond to points that have been replaced with their periodic counterpart
! and thus are not used any more
- where(rmass_inverse_elastic == 0._CUSTOM_REAL) rmass_inverse_elastic = 1._CUSTOM_REAL
- rmass_inverse_elastic(:) = 1._CUSTOM_REAL / rmass_inverse_elastic(:)
+ where(rmass_inverse_elastic_one == 0._CUSTOM_REAL) rmass_inverse_elastic_one = 1._CUSTOM_REAL
+ rmass_inverse_elastic_one(:) = 1._CUSTOM_REAL / rmass_inverse_elastic_one(:)
endif ! of if(ADD_PERIODIC_CONDITIONS)
@@ -2441,8 +2454,10 @@
allocate(veloc_elastic(3,nglob_elastic))
allocate(accel_elastic(3,nglob_elastic))
allocate(accel_elastic_adj_coupling(3,nglob_elastic))
- allocate(rmass_inverse_elastic(nglob_elastic))
+ allocate(rmass_inverse_elastic_one(nglob_elastic))
+ allocate(rmass_inverse_elastic_three(nglob_elastic))
+
if(time_stepping_scheme==2)then
allocate(displ_elastic_LDDRK(3,nglob_elastic))
allocate(veloc_elastic_LDDRK(3,nglob_elastic))
@@ -2700,7 +2715,7 @@
!---- build the global mass matrix
!
call invert_mass_matrix_init(any_elastic,any_acoustic,any_poroelastic, &
- rmass_inverse_elastic,nglob_elastic, &
+ rmass_inverse_elastic_one,nglob_elastic, &
rmass_inverse_acoustic,nglob_acoustic, &
rmass_s_inverse_poroelastic, &
rmass_w_inverse_poroelastic,nglob_poroelastic, &
@@ -2708,8 +2723,12 @@
elastic,poroelastic, &
assign_external_model,numat, &
density,poroelastcoef,porosity,tortuosity, &
- vpext,rhoext)
+ vpext,rhoext,&
+ anyabs,numabs,deltat,codeabs,rmass_inverse_elastic_three,&
+ nelemabs,vsext,xix,xiz,gammaz,gammax)
+
+
#ifdef USE_MPI
if ( nproc > 1 ) then
@@ -2774,7 +2793,7 @@
! assembling the mass matrix
call assemble_MPI_scalar(rmass_inverse_acoustic,nglob_acoustic, &
- rmass_inverse_elastic,nglob_elastic, &
+ rmass_inverse_elastic_one,nglob_elastic, &
rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic,nglob_poroelastic, &
ninterface, max_interface_size, max_ibool_interfaces_size_ac, &
max_ibool_interfaces_size_el, &
@@ -2938,11 +2957,12 @@
!---
call invert_mass_matrix(any_elastic,any_acoustic,any_poroelastic,&
- rmass_inverse_elastic,nglob_elastic, &
+ rmass_inverse_elastic_one,rmass_inverse_elastic_three,nglob_elastic, &
rmass_inverse_acoustic,nglob_acoustic, &
rmass_s_inverse_poroelastic, &
rmass_w_inverse_poroelastic,nglob_poroelastic)
+
! check the mesh, stability and number of points per wavelength
if(DISPLAY_SUBSET_OPTION == 1) then
UPPER_LIMIT_DISPLAY = nspec
@@ -3041,6 +3061,7 @@
if(j > NZ_IMAGE_color) j = NZ_IMAGE_color
iglob_image_color(i,j) = iproc
+
enddo
enddo
@@ -5107,7 +5128,7 @@
nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,mu_k,kappa_k, &
e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_rk, e11_force_rk, e13_force_rk, &
- stage_time_scheme,i_stage)
+ stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring)
if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1) then
!--- left absorbing boundary
@@ -5693,16 +5714,18 @@
if(any_elastic) then
if(time_stepping_scheme == 1)then
- accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic
- accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic
- accel_elastic(3,:) = accel_elastic(3,:) * rmass_inverse_elastic
+
+ accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic_one
+ accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic_one
+ accel_elastic(3,:) = accel_elastic(3,:) * rmass_inverse_elastic_three
+
veloc_elastic = veloc_elastic + deltatover2 * accel_elastic
endif
if(time_stepping_scheme == 2)then
- accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic
- accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic
- accel_elastic(3,:) = accel_elastic(3,:) * rmass_inverse_elastic
+ accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic_one
+ accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic_one
+ accel_elastic(3,:) = accel_elastic(3,:) * rmass_inverse_elastic_three
veloc_elastic_LDDRK = alpha_LDDRK(i_stage) * veloc_elastic_LDDRK + deltat * accel_elastic
displ_elastic_LDDRK = alpha_LDDRK(i_stage) * displ_elastic_LDDRK + deltat * veloc_elastic
@@ -5714,9 +5737,9 @@
if(time_stepping_scheme == 3)then
- accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic
- accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic
- accel_elastic(3,:) = accel_elastic(3,:) * rmass_inverse_elastic
+ accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic_one
+ accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic_one
+ accel_elastic(3,:) = accel_elastic(3,:) * rmass_inverse_elastic_three
accel_elastic_rk(1,:,i_stage) = deltat * accel_elastic(1,:)
accel_elastic_rk(2,:,i_stage) = deltat * accel_elastic(2,:)
@@ -5783,9 +5806,9 @@
endif
if(SIMULATION_TYPE == 2) then
- b_accel_elastic(1,:) = b_accel_elastic(1,:) * rmass_inverse_elastic(:)
- b_accel_elastic(2,:) = b_accel_elastic(2,:) * rmass_inverse_elastic(:)
- b_accel_elastic(3,:) = b_accel_elastic(3,:) * rmass_inverse_elastic(:)
+ b_accel_elastic(1,:) = b_accel_elastic(1,:) * rmass_inverse_elastic_one(:)
+ b_accel_elastic(2,:) = b_accel_elastic(2,:) * rmass_inverse_elastic_one(:)
+ b_accel_elastic(3,:) = b_accel_elastic(3,:) * rmass_inverse_elastic_three(:)
b_veloc_elastic = b_veloc_elastic + b_deltatover2*b_accel_elastic
endif
@@ -6631,8 +6654,8 @@
veloc_elastic(1,iglob) = veloc_elastic(1,iglob) - deltatover2*accel_elastic(1,iglob)
veloc_elastic(3,iglob) = veloc_elastic(3,iglob) - deltatover2*accel_elastic(3,iglob)
- accel_elastic(1,iglob) = accel_elastic(1,iglob) / rmass_inverse_elastic(iglob)
- accel_elastic(3,iglob) = accel_elastic(3,iglob) / rmass_inverse_elastic(iglob)
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) / rmass_inverse_elastic_one(iglob)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) / rmass_inverse_elastic_three(iglob)
! recovering original velocities and accelerations on boundaries (poro side)
velocs_poroelastic(1,iglob) = velocs_poroelastic(1,iglob) - deltatover2*accels_poroelastic(1,iglob)
velocs_poroelastic(2,iglob) = velocs_poroelastic(2,iglob) - deltatover2*accels_poroelastic(2,iglob)
@@ -6640,9 +6663,9 @@
accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) / rmass_s_inverse_poroelastic(iglob)
! assembling accelerations
accel_elastic(1,iglob) = ( accel_elastic(1,iglob) + accels_poroelastic(1,iglob) ) / &
- ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+ ( 1.0/rmass_inverse_elastic_one(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
accel_elastic(3,iglob) = ( accel_elastic(3,iglob) + accels_poroelastic(2,iglob) ) / &
- ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+ ( 1.0/rmass_inverse_elastic_three(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
accels_poroelastic(1,iglob) = accel_elastic(1,iglob)
accels_poroelastic(2,iglob) = accel_elastic(3,iglob)
! updating velocities
@@ -6880,8 +6903,8 @@
if(SIMULATION_TYPE == 2) then
b_veloc_elastic(1,iglob) = b_veloc_elastic(1,iglob) - b_deltatover2*b_accel_elastic(1,iglob)
b_veloc_elastic(3,iglob) = b_veloc_elastic(3,iglob) - b_deltatover2*b_accel_elastic(3,iglob)
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) / rmass_inverse_elastic(iglob)
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) / rmass_inverse_elastic(iglob)
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) / rmass_inverse_elastic_one(iglob)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) / rmass_inverse_elastic_three(iglob)
! recovering original velocities and accelerations on boundaries (poro side)
b_velocs_poroelastic(1,iglob) = b_velocs_poroelastic(1,iglob) - b_deltatover2*b_accels_poroelastic(1,iglob)
b_velocs_poroelastic(2,iglob) = b_velocs_poroelastic(2,iglob) - b_deltatover2*b_accels_poroelastic(2,iglob)
@@ -6889,9 +6912,9 @@
b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) / rmass_s_inverse_poroelastic(iglob)
! assembling accelerations
b_accel_elastic(1,iglob) = ( b_accel_elastic(1,iglob) + b_accels_poroelastic(1,iglob) ) / &
- ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+ ( 1.0/rmass_inverse_elastic_one(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
b_accel_elastic(3,iglob) = ( b_accel_elastic(3,iglob) + b_accels_poroelastic(2,iglob) ) / &
- ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+ ( 1.0/rmass_inverse_elastic_three(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
b_accels_poroelastic(1,iglob) = b_accel_elastic(1,iglob)
b_accels_poroelastic(2,iglob) = b_accel_elastic(3,iglob)
! updating velocities
More information about the CIG-COMMITS
mailing list