diff --git a/DATA/Par_file b/DATA/Par_file index 3119e589e..a989697dd 100644 --- a/DATA/Par_file +++ b/DATA/Par_file @@ -268,6 +268,10 @@ USE_FORCE_POINT_SOURCE = .false. # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +# uses a sin-squared STF useful for PEGS calculations +# see https://doi.org/10.1016/j.epsl.2020.116150 for definition +USE_SINSQ_STF = .false. + # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/benchmarks/attenuation_benchmark_GJI_2002_versus_normal_modes/DATA/Par_file b/EXAMPLES/benchmarks/attenuation_benchmark_GJI_2002_versus_normal_modes/DATA/Par_file index ef3b2bc3b..fbdaad588 100644 --- a/EXAMPLES/benchmarks/attenuation_benchmark_GJI_2002_versus_normal_modes/DATA/Par_file +++ b/EXAMPLES/benchmarks/attenuation_benchmark_GJI_2002_versus_normal_modes/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/benchmarks/prem_ani_noocean_Q_Cowling_Bolivia_9s/SEMD/DATA/Par_file b/EXAMPLES/benchmarks/prem_ani_noocean_Q_Cowling_Bolivia_9s/SEMD/DATA/Par_file index b261874a8..591123211 100644 --- a/EXAMPLES/benchmarks/prem_ani_noocean_Q_Cowling_Bolivia_9s/SEMD/DATA/Par_file +++ b/EXAMPLES/benchmarks/prem_ani_noocean_Q_Cowling_Bolivia_9s/SEMD/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/benchmarks/prem_ani_noocean_Q_Cowling_Vanuatu/SEMD/DATA/Par_file b/EXAMPLES/benchmarks/prem_ani_noocean_Q_Cowling_Vanuatu/SEMD/DATA/Par_file index 84c8a3b88..93316829a 100644 --- a/EXAMPLES/benchmarks/prem_ani_noocean_Q_Cowling_Vanuatu/SEMD/DATA/Par_file +++ b/EXAMPLES/benchmarks/prem_ani_noocean_Q_Cowling_Vanuatu/SEMD/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/benchmarks/prem_iso_noocean_Q_Cowling_Vanuatu/SEMD/DATA/Par_file b/EXAMPLES/benchmarks/prem_iso_noocean_Q_Cowling_Vanuatu/SEMD/DATA/Par_file index 11de9bd1a..7d7add899 100644 --- a/EXAMPLES/benchmarks/prem_iso_noocean_Q_Cowling_Vanuatu/SEMD/DATA/Par_file +++ b/EXAMPLES/benchmarks/prem_iso_noocean_Q_Cowling_Vanuatu/SEMD/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/benchmarks/undo_attenuation_benchmark_GJI_2016/DATA/Par_file b/EXAMPLES/benchmarks/undo_attenuation_benchmark_GJI_2016/DATA/Par_file index 4844a60fe..ee0bcd5be 100644 --- a/EXAMPLES/benchmarks/undo_attenuation_benchmark_GJI_2016/DATA/Par_file +++ b/EXAMPLES/benchmarks/undo_attenuation_benchmark_GJI_2016/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/global_PREM_kernels/amplitude/DATA/Par_file b/EXAMPLES/global_PREM_kernels/amplitude/DATA/Par_file index 859b43e29..a5571b003 100644 --- a/EXAMPLES/global_PREM_kernels/amplitude/DATA/Par_file +++ b/EXAMPLES/global_PREM_kernels/amplitude/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/global_PREM_kernels/traveltime/DATA/Par_file b/EXAMPLES/global_PREM_kernels/traveltime/DATA/Par_file index 859b43e29..a5571b003 100644 --- a/EXAMPLES/global_PREM_kernels/traveltime/DATA/Par_file +++ b/EXAMPLES/global_PREM_kernels/traveltime/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/global_s362ani_shakemovie/DATA/Par_file b/EXAMPLES/global_s362ani_shakemovie/DATA/Par_file index 83774589d..4e4d240e4 100644 --- a/EXAMPLES/global_s362ani_shakemovie/DATA/Par_file +++ b/EXAMPLES/global_s362ani_shakemovie/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/global_small/DATA/Par_file b/EXAMPLES/global_small/DATA/Par_file index 5982ad6b6..d29bcc5c6 100644 --- a/EXAMPLES/global_small/DATA/Par_file +++ b/EXAMPLES/global_small/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/mars_global/DATA/Par_file b/EXAMPLES/mars_global/DATA/Par_file index 36ac4cdd2..45328e4aa 100644 --- a/EXAMPLES/mars_global/DATA/Par_file +++ b/EXAMPLES/mars_global/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/mars_regional/DATA/Par_file b/EXAMPLES/mars_regional/DATA/Par_file index f075b8939..c3a607f0e 100644 --- a/EXAMPLES/mars_regional/DATA/Par_file +++ b/EXAMPLES/mars_regional/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/moon_global/DATA/Par_file b/EXAMPLES/moon_global/DATA/Par_file index a8b9e183f..434fd8e64 100644 --- a/EXAMPLES/moon_global/DATA/Par_file +++ b/EXAMPLES/moon_global/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_1_attenuation b/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_1_attenuation index 04078ad8e..9d54381d7 100644 --- a/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_1_attenuation +++ b/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_1_attenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_1_noattenuation b/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_1_noattenuation index ba5ac5bb4..375b0738b 100644 --- a/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_1_noattenuation +++ b/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_1_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_2_attenuation b/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_2_attenuation index ef8ff8b29..d4329a613 100644 --- a/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_2_attenuation +++ b/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_2_attenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_2_noattenuation b/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_2_noattenuation index 341952fa1..1a9d918d2 100644 --- a/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_2_noattenuation +++ b/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_2_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_3_noattenuation b/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_3_noattenuation index bffd5329e..4892c1777 100644 --- a/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_3_noattenuation +++ b/EXAMPLES/noise_examples/global_long/DATA/Par_file_NOISE_3_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_1_attenuation b/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_1_attenuation index 6ad0bb5d1..b8e5961e5 100644 --- a/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_1_attenuation +++ b/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_1_attenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_1_noattenuation b/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_1_noattenuation index 324b10fff..85e6dab49 100644 --- a/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_1_noattenuation +++ b/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_1_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_2_attenuation b/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_2_attenuation index c5554e2d6..56b8566ad 100644 --- a/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_2_attenuation +++ b/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_2_attenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_2_noattenuation b/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_2_noattenuation index a5ae133ac..eccf0a1a9 100644 --- a/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_2_noattenuation +++ b/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_2_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_3_noattenuation b/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_3_noattenuation index 764d2609b..3bcb8695e 100644 --- a/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_3_noattenuation +++ b/EXAMPLES/noise_examples/global_short/DATA/Par_file_NOISE_3_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_1_attenuation b/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_1_attenuation index ae4eed042..5a402dbf8 100644 --- a/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_1_attenuation +++ b/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_1_attenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_1_noattenuation b/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_1_noattenuation index a378e79fa..2e2e9da33 100644 --- a/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_1_noattenuation +++ b/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_1_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_2_attenuation b/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_2_attenuation index dee8b2c07..033ca4c20 100644 --- a/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_2_attenuation +++ b/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_2_attenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_2_noattenuation b/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_2_noattenuation index d233830dd..5ec06719a 100644 --- a/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_2_noattenuation +++ b/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_2_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_3_noattenuation b/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_3_noattenuation index 70b326358..23561ea03 100644 --- a/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_3_noattenuation +++ b/EXAMPLES/noise_examples/regional/DATA/Par_file_NOISE_3_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_1_attenuation b/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_1_attenuation index e83f80e1c..2601048c3 100644 --- a/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_1_attenuation +++ b/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_1_attenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_1_noattenuation b/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_1_noattenuation index 00c27b79b..8c037bcd6 100644 --- a/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_1_noattenuation +++ b/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_1_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_2_attenuation b/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_2_attenuation index 38f346793..b328eb40f 100644 --- a/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_2_attenuation +++ b/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_2_attenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_2_noattenuation b/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_2_noattenuation index 98286ebcd..d0db3d94d 100644 --- a/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_2_noattenuation +++ b/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_2_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_3_noattenuation b/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_3_noattenuation index 266b3bc5d..81d00e991 100644 --- a/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_3_noattenuation +++ b/EXAMPLES/noise_examples/test_global/DATA/Par_file_NOISE_3_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_1_attenuation b/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_1_attenuation index e874fc23d..7e0c3f49b 100644 --- a/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_1_attenuation +++ b/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_1_attenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_1_noattenuation b/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_1_noattenuation index 5d5c8b758..8d8f68edd 100644 --- a/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_1_noattenuation +++ b/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_1_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_2_attenuation b/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_2_attenuation index 824fa401a..16733cbff 100644 --- a/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_2_attenuation +++ b/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_2_attenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_2_noattenuation b/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_2_noattenuation index 0f16714af..a5de0a029 100644 --- a/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_2_noattenuation +++ b/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_2_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_3_noattenuation b/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_3_noattenuation index 0cb7ae708..89366ef9c 100644 --- a/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_3_noattenuation +++ b/EXAMPLES/noise_examples/test_regional/DATA/Par_file_NOISE_3_noattenuation @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/point_force/DATA/Par_file b/EXAMPLES/point_force/DATA/Par_file index b9c7cebcf..ca70ab4b8 100644 --- a/EXAMPLES/point_force/DATA/Par_file +++ b/EXAMPLES/point_force/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .true. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/regional_Berkeley/DATA/Par_file b/EXAMPLES/regional_Berkeley/DATA/Par_file index 7415d4edf..8afa9025a 100644 --- a/EXAMPLES/regional_Berkeley/DATA/Par_file +++ b/EXAMPLES/regional_Berkeley/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/regional_Berkeley/REF_SEIS/Par_file b/EXAMPLES/regional_Berkeley/REF_SEIS/Par_file index 7415d4edf..8afa9025a 100644 --- a/EXAMPLES/regional_Berkeley/REF_SEIS/Par_file +++ b/EXAMPLES/regional_Berkeley/REF_SEIS/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/regional_EMC_model/DATA/Par_file b/EXAMPLES/regional_EMC_model/DATA/Par_file index 7665801a5..a003af120 100644 --- a/EXAMPLES/regional_EMC_model/DATA/Par_file +++ b/EXAMPLES/regional_EMC_model/DATA/Par_file @@ -287,6 +287,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/regional_EMC_model/DATA/Par_file.SCEC-CVM b/EXAMPLES/regional_EMC_model/DATA/Par_file.SCEC-CVM index 008d9d3fe..d540d0581 100644 --- a/EXAMPLES/regional_EMC_model/DATA/Par_file.SCEC-CVM +++ b/EXAMPLES/regional_EMC_model/DATA/Par_file.SCEC-CVM @@ -287,6 +287,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/regional_Greece_noise_small/DATA/Par_file b/EXAMPLES/regional_Greece_noise_small/DATA/Par_file index 12e098c35..77f5ff436 100644 --- a/EXAMPLES/regional_Greece_noise_small/DATA/Par_file +++ b/EXAMPLES/regional_Greece_noise_small/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/regional_Greece_small/DATA/Par_file b/EXAMPLES/regional_Greece_small/DATA/Par_file index 1ce8c7f39..d01c69453 100644 --- a/EXAMPLES/regional_Greece_small/DATA/Par_file +++ b/EXAMPLES/regional_Greece_small/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/regional_Greece_small/OUTPUT_FILES/Par_file b/EXAMPLES/regional_Greece_small/OUTPUT_FILES/Par_file new file mode 100644 index 000000000..242fdf449 --- /dev/null +++ b/EXAMPLES/regional_Greece_small/OUTPUT_FILES/Par_file @@ -0,0 +1,440 @@ +#----------------------------------------------------------- +# +# Simulation input parameters +# +#----------------------------------------------------------- + +# forward or adjoint simulation +SIMULATION_TYPE = 1 # set to 1 for forward simulations, 2 for adjoint simulations for sources, and 3 for kernel simulations +NOISE_TOMOGRAPHY = 0 # flag of noise tomography, three steps (1,2,3). If earthquake simulation, set it to 0. +SAVE_FORWARD = .false. # save last frame of forward simulation or not + +# number of chunks (1,2,3 or 6) +NCHUNKS = 1 + +# angular width of the first chunk (not used if full sphere with six chunks) +ANGULAR_WIDTH_XI_IN_DEGREES = 20.d0 # angular size of a chunk +ANGULAR_WIDTH_ETA_IN_DEGREES = 20.d0 +CENTER_LATITUDE_IN_DEGREES = 40.d0 +CENTER_LONGITUDE_IN_DEGREES = 25.d0 +GAMMA_ROTATION_AZIMUTH = 0.d0 + +# number of elements at the surface along the two sides of the first chunk +# (must be multiple of 16 and 8 * multiple of NPROC below) +NEX_XI = 64 +NEX_ETA = 64 + +# number of MPI processors along the two sides of the first chunk +NPROC_XI = 2 +NPROC_ETA = 2 + +#----------------------------------------------------------- +# +# Model +# +#----------------------------------------------------------- + +# 1D models with real structure: +# 1D_isotropic_prem, 1D_transversely_isotropic_prem, 1D_iasp91, 1D_1066a, 1D_ak135f_no_mud, 1D_ref, 1D_ref_iso, 1D_jp3d, 1D_sea99, 1D_Berkeley +# +# 1D models with only one fictitious averaged crustal layer: +# 1D_isotropic_prem_onecrust, 1D_transversely_isotropic_prem_onecrust, 1D_iasp91_onecrust, 1D_1066a_onecrust, 1D_ak135f_no_mud_onecrust +# +# fully 3D models: +# transversely_isotropic_prem_plus_3D_crust_2.0, 3D_anisotropic, 3D_attenuation, +# s20rts, s40rts, s362ani, s362iso, s362wmani, s362ani_prem, s362ani_3DQ, s362iso_3DQ, +# s29ea, sea99_jp3d1994, sea99, jp3d1994, heterogen, full_sh, sgloberani_aniso, sgloberani_iso, +# spiral, emc_model, emc_model_qmu, emc_model_tiso, emc_model_tiso_qmu, semucb_A3d, semucb_A3d_3dQ +# (see manual for more...) +# +# 3D crustal models: +# crust1.0, crust2.0, EPcrust, EuCRUST, crustmaps, crustSH +# +# Mars models: +# 1D_Sohl, 1D_Sohl_3D_crust, 1D_case65TAY, 1D_case65TAY_3D_crust, mars_1D, mars_1D_3D_crust +# +# Moon models: +# vpremoon +# +# 3D models with 3D crust: append "_**crustname**" to the mantle model name +# to take a 3D crustal model (by default crust2.0 is taken for 3D mantle models) +# e.g. s20rts_crust1.0, s362ani_crustmaps, full_sh_crustSH, sglobe_EPcrust, etc. +# +# 3D models with 1D crust: append "_1Dcrust" to the 3D model name +# to take the 1D crustal model from the +# associated reference model rather than the default 3D crustal model +# e.g. s20rts_1Dcrust, s362ani_1Dcrust, etc. +# +MODEL = 1D_transversely_isotropic_prem + +# parameters describing the Earth model +OCEANS = .true. +ELLIPTICITY = .true. +TOPOGRAPHY = .true. +GRAVITY = .true. +ROTATION = .true. +ATTENUATION = .true. + +# full gravity calculation by solving Poisson's equation for gravity potential instead of using a Cowling approximation +# (must have also GRAVITY flag set to .true. to become active) +FULL_GRAVITY = .false. +# for full gravity calculation, set to 0 == builtin or 1 == PETSc Poisson solver +# (the PETSc solver option needs the PETSc library installed; code configuration --with-petsc) +POISSON_SOLVER = 0 + +# record length in minutes +RECORD_LENGTH_IN_MINUTES = 2.5d0 + +#----------------------------------------------------------- +# +# Mesh +# +#----------------------------------------------------------- + +## regional mesh cut-off +# using this flag will cut-off the mesh in the mantle at a layer matching to the given cut-off depth. +# this flag only has an effect for regional simulations, i.e., for NCHUNKS values less than 6. +REGIONAL_MESH_CUTOFF = .false. + +# regional mesh cut-off depth (in km) +# possible selections are: 24.4d0, 80.d0, 220.d0, 400.d0, 600.d0, 670.d0, 771.d0 +REGIONAL_MESH_CUTOFF_DEPTH = 400.d0 + +# regional mesh cut-off w/ a second doubling layer below 220km interface +# (by default, a first doubling layer will be added below the Moho, and a second one below the 771km-depth layer. +# Setting this flag to .true., will move the second one below the 220km-depth layer for regional mesh cut-offs only.) +REGIONAL_MESH_ADD_2ND_DOUBLING = .false. + +#----------------------------------------------------------- +# +# Absorbing boundary conditions +# +#----------------------------------------------------------- + +# absorbing boundary conditions for a regional simulation +ABSORBING_CONDITIONS = .true. + +# run global simulation for a circular region and apply high attenuation for the rest of the model +# this creates an absorbing boundary with less reflection than stacey at a cost of 6x the computational cost +# NCHUNKS must be set to 6 to enable this flag +ABSORB_USING_GLOBAL_SPONGE = .false. + +# location and size of the region with no sponge (the region to run simulation) +SPONGE_LATITUDE_IN_DEGREES = 40.d0 +SPONGE_LONGITUDE_IN_DEGREES = 25.d0 +SPONGE_RADIUS_IN_DEGREES = 25.d0 + +#----------------------------------------------------------- +# +# undoing attenuation for sensitivity kernel calculations +# +#----------------------------------------------------------- + +# to undo attenuation for sensitivity kernel calculations or forward runs with SAVE_FORWARD +# use one (and only one) of the two flags below. UNDO_ATTENUATION is much better (it is exact) +# but requires a significant amount of disk space for temporary storage. +PARTIAL_PHYS_DISPERSION_ONLY = .true. +UNDO_ATTENUATION = .false. + +## undo attenuation memory +# How much memory (in GB) is installed on your machine per CPU core +# (only used for UNDO_ATTENUATION, can be ignored otherwise) +# Beware, this value MUST be given per core, i.e. per MPI thread, i.e. per MPI rank, NOT per node. +# This value is for instance: +# - 4 GB on Tiger at Princeton +# - 4 GB on TGCC Curie in Paris +# - 4 GB on Titan at ORNL when using CPUs only (no GPUs); start your run with "aprun -n$NPROC -N8 -S4 -j1" +# - 2 GB on the machine used by Christina Morency +# - 2 GB on the TACC machine used by Min Chen +# - 1.5 GB on the GPU cluster in Marseille +# When running on GPU machines, it is simpler to set PERCENT_OF_MEM_TO_USE_PER_CORE = 100.d0 +# and then set MEMORY_INSTALLED_PER_CORE_IN_GB to the amount of memory that you estimate is free (rather than installed) +# on the host of the GPU card while running your GPU job. +# For GPU runs on Titan at ORNL, use PERCENT_OF_MEM_TO_USE_PER_CORE = 100.d0 and MEMORY_INSTALLED_PER_CORE_IN_GB = 25.d0 +# and run your job with "aprun -n$NPROC -N1 -S1 -j1" +# (each host has 32 GB on Titan, each GPU has 6 GB, thus even if all the GPU arrays are duplicated on the host +# this leaves 32 - 6 = 26 GB free on the host; leaving 1 GB for the Linux system, we can safely use 100% of 25 GB) +MEMORY_INSTALLED_PER_CORE_IN_GB = 4.d0 +# What percentage of this total do you allow us to use for arrays to undo attenuation, keeping in mind that you +# need to leave some memory available for the GNU/Linux system to run +# (a typical value is 85%; any value below is fine but the code will then save a lot of data to disk; +# values above, say 90% or 92%, can be OK on some systems but can make the adjoint code run out of memory +# on other systems, depending on how much memory per node the GNU/Linux system needs for itself; thus you can try +# a higher value and if the adjoint crashes then try again with a lower value) +PERCENT_OF_MEM_TO_USE_PER_CORE = 85.d0 + +## exact mass matrices for rotation +# three mass matrices instead of one are needed to handle rotation very accurately; +# otherwise rotation is handled slightly less accurately (but still reasonably well); +# set to .true. if you are interested in precise effects related to rotation; +# set to .false. if you are solving very large inverse problems at high frequency and also undoing attenuation exactly +# +# using the UNDO_ATTENUATION flag above, in which case saving as much memory as possible can be a good idea. +# You can also safely set it to .false. if you are not in a period range in which rotation matters, +# e.g. if you are targetting very short-period body waves. if in doubt, set to .true. +# +# You can safeely set it to .true. if you have ABSORBING_CONDITIONS above, because in that case the code +# will use three mass matrices anyway and thus there is no additional memory cost. +# this flag is of course unused if ROTATION above is set to .false. +EXACT_MASS_MATRIX_FOR_ROTATION = .false. + +#----------------------------------------------------------- +# +# LDDRK time scheme +# +#----------------------------------------------------------- + +# this for LDDRK high-order time scheme instead of Newmark +USE_LDDRK = .false. + +# the maximum CFL of LDDRK is significantly higher than that of the Newmark scheme, +# in a ratio that is theoretically 1.327 / 0.697 = 1.15 / 0.604 = 1.903 for a solid with Poisson's ratio = 0.25 +# and for a fluid (see the manual of the 2D code, SPECFEM2D, Tables 4.1 and 4.2, and that ratio does not +# depend on whether we are in 2D or in 3D). However in practice a ratio of about 1.5 to 1.7 is often safer +# (for instance for models with a large range of Poisson's ratio values). +# Since the code computes the time step using the Newmark scheme, for LDDRK we will simply +# multiply that time step by this ratio when LDDRK is on and when flag INCREASE_CFL_FOR_LDDRK is true. +INCREASE_CFL_FOR_LDDRK = .true. +RATIO_BY_WHICH_TO_INCREASE_IT = 1.5d0 + +#----------------------------------------------------------- +# +# Visualization +# +#----------------------------------------------------------- + +# save AVS or OpenDX movies +#MOVIE_COARSE saves movie only at corners of elements (SURFACE OR VOLUME) +#MOVIE_COARSE does not work with create_movie_AVS_DX +MOVIE_SURFACE = .false. +MOVIE_VOLUME = .false. +MOVIE_COARSE = .true. +NTSTEP_BETWEEN_FRAMES = 50 +HDUR_MOVIE = 0.d0 + +# save movie in volume. Will save element if center of element is in prescribed volume +# top/bottom: depth in KM, use MOVIE_TOP = -100 to make sure the surface is stored. +# west/east: longitude, degrees East [-180/180] top/bottom: latitute, degrees North [-90/90] +# start/stop: frames will be stored at MOVIE_START + i*NSTEP_BETWEEN_FRAMES, where i=(0,1,2..) and iNSTEP_BETWEEN_FRAMES <= MOVIE_STOP +# movie_volume_type: 1=strain, 2=time integral of strain, 3=\mu*time integral of strain +# type 4 saves the trace and deviatoric stress in the whole volume, 5=displacement, 6=velocity +MOVIE_VOLUME_TYPE = 2 +MOVIE_TOP_KM = -100.0 +MOVIE_BOTTOM_KM = 1000.0 +MOVIE_WEST_DEG = -90.0 +MOVIE_EAST_DEG = 90.0 +MOVIE_NORTH_DEG = 90.0 +MOVIE_SOUTH_DEG = -90.0 +MOVIE_START = 0 +MOVIE_STOP = 40000 + +# save mesh files to check the mesh +SAVE_MESH_FILES = .true. + +# restart files (number of runs can be 1 or higher, choose 1 for no restart files) +NUMBER_OF_RUNS = 1 +NUMBER_OF_THIS_RUN = 1 + +# path to store the local database files on each node +LOCAL_PATH = ./DATABASES_MPI +# temporary wavefield/kernel/movie files +LOCAL_TMP_PATH = ./DATABASES_MPI + +# interval at which we output time step info and max of norm of displacement +NTSTEP_BETWEEN_OUTPUT_INFO = 200 + +#----------------------------------------------------------- +# +# Sources +# +#----------------------------------------------------------- + +# use a (tilted) FORCESOLUTION force point source (or several) instead of a CMTSOLUTION moment-tensor source. +# This can be useful e.g. for asteroid simulations +# in which the source is a vertical force, normal force, tilted force, impact etc. +# If this flag is turned on, the FORCESOLUTION file must be edited by giving: +# - the corresponding time-shift parameter, +# - the half duration parameter of the source, +# - the coordinates of the source, +# - the source time function of the source, +# - the magnitude of the force source, +# - the components of a (non necessarily unitary) direction vector for the force source in the E/N/Z_UP basis. +# The direction vector is made unitary internally in the code and thus only its direction matters here; +# its norm is ignored and the norm of the force used is the factor force source times the source time function. +USE_FORCE_POINT_SOURCE = .false. + +# use monochromatic source time function for CMTSOLUTION moment-tensor source. +# half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format +# default is .false. which uses Heaviside function +USE_MONOCHROMATIC_CMT_SOURCE = .false. + + +# print source time function +PRINT_SOURCE_TIME_FUNCTION = .true. + +#----------------------------------------------------------- +# +# Seismograms +# +#----------------------------------------------------------- + +# interval in time steps for temporary writing of seismograms +NTSTEP_BETWEEN_OUTPUT_SEISMOS = 5000000 + +# set to n to reduce the sampling rate of output seismograms by a factor of n +# defaults to 1, which means no down-sampling +NTSTEP_BETWEEN_OUTPUT_SAMPLE = 1 + +# option to save strain seismograms +# this option is useful for strain Green's tensor +# this feature is currently under development +SAVE_SEISMOGRAMS_STRAIN = .false. + +# save seismograms also when running the adjoint runs for an inverse problem +# (usually they are unused and not very meaningful, leave this off in almost all cases) +SAVE_SEISMOGRAMS_IN_ADJOINT_RUN = .true. + +# output format for the seismograms (one can use either or all of the three formats) +OUTPUT_SEISMOS_ASCII_TEXT = .true. +OUTPUT_SEISMOS_SAC_ALPHANUM = .false. +OUTPUT_SEISMOS_SAC_BINARY = .false. +OUTPUT_SEISMOS_ASDF = .false. +OUTPUT_SEISMOS_3D_ARRAY = .false. +OUTPUT_SEISMOS_HDF5 = .false. + +# rotate seismograms to Radial-Transverse-Z or use default North-East-Z reference frame +ROTATE_SEISMOGRAMS_RT = .false. + +# decide if main process writes all the seismograms or if all processes do it in parallel +WRITE_SEISMOGRAMS_BY_MAIN = .false. + +# save all seismograms in one large combined file instead of one file per seismogram +# to avoid overloading shared non-local file systems such as LUSTRE or GPFS for instance +SAVE_ALL_SEISMOS_IN_ONE_FILE = .false. +USE_BINARY_FOR_LARGE_FILE = .false. + +# flag to impose receivers at the surface or allow them to be buried +RECEIVERS_CAN_BE_BURIED = .true. + +#----------------------------------------------------------- +# +# Adjoint kernel outputs +# +#----------------------------------------------------------- + +# interval in time steps for reading adjoint traces +# 0 = read the whole adjoint sources at start time +NTSTEP_BETWEEN_READ_ADJSRC = 10000 + +# use ASDF format for reading the adjoint sources +READ_ADJSRC_ASDF = .false. + +# this parameter must be set to .true. to compute anisotropic kernels +# in crust and mantle (related to the 21 Cij in geographical coordinates) +# default is .false. to compute isotropic kernels (related to alpha and beta) +ANISOTROPIC_KL = .false. + +# output only transverse isotropic kernels (alpha_v,alpha_h,beta_v,beta_h,eta,rho) +# rather than fully anisotropic kernels when ANISOTROPIC_KL above is set to .true. +# means to save radial anisotropic kernels, i.e., sensitivity kernels for beta_v, beta_h, etc. +SAVE_TRANSVERSE_KL_ONLY = .false. + +# output only the kernels used for the current azimuthally anisotropic inversions of surface waves, +# i.e., bulk_c, bulk_betav, bulk_betah, eta, Gc_prime, Gs_prime and rho +# (Gc' & Gs' which are the normalized Gc & Gs kernels by isotropic \rho\beta of the 1D reference model) +SAVE_AZIMUTHAL_ANISO_KL_ONLY = .false. + +# output approximate Hessian in crust mantle region. +# means to save the preconditioning for gradients, they are cross correlations between forward and adjoint accelerations. +APPROXIMATE_HESS_KL = .false. + +# forces transverse isotropy for all mantle elements +# (default is to use transverse isotropy only between crust and 220) +# means we allow radial anisotropy throughout the whole crust/mantle region +USE_FULL_TISO_MANTLE = .false. + +# output kernel mask to zero out source region +# to remove large values near the sources in the sensitivity kernels +SAVE_SOURCE_MASK = .false. + +# output kernels on a regular grid instead of on the GLL mesh points (a bit expensive) +SAVE_REGULAR_KL = .false. + +# compute steady state kernels for source encoding +STEADY_STATE_KERNEL = .false. +STEADY_STATE_LENGTH_IN_MINUTES = 0.d0 + +#----------------------------------------------------------- + +# Dimitri Komatitsch, July 2014, CNRS Marseille, France: +# added the ability to run several calculations (several earthquakes) +# in an embarrassingly-parallel fashion from within the same run; +# this can be useful when using a very large supercomputer to compute +# many earthquakes in a catalog, in which case it can be better from +# a batch job submission point of view to start fewer and much larger jobs, +# each of them computing several earthquakes in parallel. +# To turn that option on, set parameter NUMBER_OF_SIMULTANEOUS_RUNS to a value greater than 1. +# To implement that, we create NUMBER_OF_SIMULTANEOUS_RUNS MPI sub-communicators, +# each of them being labeled "my_local_mpi_comm_world", and we use them +# in all the routines in "src/shared/parallel.f90", except in MPI_ABORT() because in that case +# we need to kill the entire run. +# When that option is on, of course the number of processor cores used to start +# the code in the batch system must be a multiple of NUMBER_OF_SIMULTANEOUS_RUNS, +# all the individual runs must use the same number of processor cores, +# which as usual is NPROC in the Par_file, +# and thus the total number of processor cores to request from the batch system +# should be NUMBER_OF_SIMULTANEOUS_RUNS * NPROC. +# All the runs to perform must be placed in directories called run0001, run0002, run0003 and so on +# (with exactly four digits). +# +# Imagine you have 10 independent calculations to do, each of them on 100 cores; you have three options: +# +# 1/ submit 10 jobs to the batch system +# +# 2/ submit a single job on 1000 cores to the batch, and in that script create a sub-array of jobs to start 10 jobs, +# each running on 100 cores (see e.g. http://www.schedmd.com/slurmdocs/job_array.html ) +# +# 3/ submit a single job on 1000 cores to the batch, start SPECFEM3D on 1000 cores, create 10 sub-communicators, +# cd into one of 10 subdirectories (called e.g. run0001, run0002,... run0010) depending on the sub-communicator +# your MPI rank belongs to, and run normally on 100 cores using that sub-communicator. +# +# The option below implements 3/. +# +NUMBER_OF_SIMULTANEOUS_RUNS = 1 + +# if we perform simultaneous runs in parallel, if only the source and receivers vary between these runs +# but not the mesh nor the model (velocity and density) then we can also read the mesh and model files +# from a single run in the beginning and broadcast them to all the others; for a large number of simultaneous +# runs for instance when solving inverse problems iteratively this can DRASTICALLY reduce I/Os to disk in the solver +# (by a factor equal to NUMBER_OF_SIMULTANEOUS_RUNS), and reducing I/Os is crucial in the case of huge runs. +# Thus, always set this option to .true. if the mesh and the model are the same for all simultaneous runs. +# In that case there is no need to duplicate the mesh and model file database (the content of the DATABASES_MPI +# directories) in each of the run0001, run0002,... directories, it is sufficient to have one in run0001 +# and the code will broadcast it to the others) +BROADCAST_SAME_MESH_AND_MODEL = .false. + +#----------------------------------------------------------- + +# set to true to use GPUs +GPU_MODE = .false. +# Only used if GPU_MODE = .true. : +GPU_RUNTIME = 1 +# 2 (OpenCL), 1 (Cuda) ou 0 (Compile-time -- does not work if configured with --with-cuda *AND* --with-opencl) +GPU_PLATFORM = NVIDIA +GPU_DEVICE = Tesla + +# set to true to use the ADIOS library for I/Os +ADIOS_ENABLED = .false. +ADIOS_FOR_FORWARD_ARRAYS = .true. +ADIOS_FOR_MPI_ARRAYS = .true. +ADIOS_FOR_ARRAYS_SOLVER = .true. +ADIOS_FOR_SOLVER_MESHFILES = .true. +ADIOS_FOR_AVS_DX = .true. +ADIOS_FOR_KERNELS = .true. +ADIOS_FOR_MODELS = .true. +ADIOS_FOR_UNDO_ATTENUATION = .true. + +# HDF5 Database I/O +# (note the flags for HDF5 and ADIOS are mutually exclusive, only one can be used) +HDF5_ENABLED = .false. + diff --git a/EXAMPLES/regional_Greece_small_LDDRK/DATA/Par_file b/EXAMPLES/regional_Greece_small_LDDRK/DATA/Par_file index fdbd2104a..8c33013d2 100644 --- a/EXAMPLES/regional_Greece_small_LDDRK/DATA/Par_file +++ b/EXAMPLES/regional_Greece_small_LDDRK/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/regional_Mexico_noise_non_uniform/DATA/Par_file b/EXAMPLES/regional_Mexico_noise_non_uniform/DATA/Par_file index 2054c3e27..b1671462f 100644 --- a/EXAMPLES/regional_Mexico_noise_non_uniform/DATA/Par_file +++ b/EXAMPLES/regional_Mexico_noise_non_uniform/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/regional_Mexico_noise_non_uniform/DATA/Par_file_step1 b/EXAMPLES/regional_Mexico_noise_non_uniform/DATA/Par_file_step1 index 67851ff63..3d219be60 100644 --- a/EXAMPLES/regional_Mexico_noise_non_uniform/DATA/Par_file_step1 +++ b/EXAMPLES/regional_Mexico_noise_non_uniform/DATA/Par_file_step1 @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/regional_Mexico_noise_non_uniform/DATA/Par_file_step2 b/EXAMPLES/regional_Mexico_noise_non_uniform/DATA/Par_file_step2 index 2054c3e27..b1671462f 100644 --- a/EXAMPLES/regional_Mexico_noise_non_uniform/DATA/Par_file_step2 +++ b/EXAMPLES/regional_Mexico_noise_non_uniform/DATA/Par_file_step2 @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/regional_MiddleEast/DATA/Par_file b/EXAMPLES/regional_MiddleEast/DATA/Par_file index e0581e1e1..3461b266a 100644 --- a/EXAMPLES/regional_MiddleEast/DATA/Par_file +++ b/EXAMPLES/regional_MiddleEast/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/regional_s40rts/DATA/Par_file b/EXAMPLES/regional_s40rts/DATA/Par_file index 07a18ea7b..bb922c58b 100644 --- a/EXAMPLES/regional_s40rts/DATA/Par_file +++ b/EXAMPLES/regional_s40rts/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/regional_sgloberani/DATA/Par_file b/EXAMPLES/regional_sgloberani/DATA/Par_file index 0e5629023..5b043a3fd 100644 --- a/EXAMPLES/regional_sgloberani/DATA/Par_file +++ b/EXAMPLES/regional_sgloberani/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/regional_simultaneous_runs/DATA/Par_file b/EXAMPLES/regional_simultaneous_runs/DATA/Par_file index c8bf51e2b..710ce94a8 100644 --- a/EXAMPLES/regional_simultaneous_runs/DATA/Par_file +++ b/EXAMPLES/regional_simultaneous_runs/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/regional_small_adjoint/DATA/Par_file b/EXAMPLES/regional_small_adjoint/DATA/Par_file index 8271b8a85..43350b7a5 100644 --- a/EXAMPLES/regional_small_adjoint/DATA/Par_file +++ b/EXAMPLES/regional_small_adjoint/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/regular_kernel/DATA/Par_file b/EXAMPLES/regular_kernel/DATA/Par_file index b0472d4a0..29186a1d5 100644 --- a/EXAMPLES/regular_kernel/DATA/Par_file +++ b/EXAMPLES/regular_kernel/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. diff --git a/EXAMPLES/small_benchmark_run_to_test_more_complex_Earth/DATA/Par_file b/EXAMPLES/small_benchmark_run_to_test_more_complex_Earth/DATA/Par_file index 0afcb81d5..ab3fb9b0b 100644 --- a/EXAMPLES/small_benchmark_run_to_test_more_complex_Earth/DATA/Par_file +++ b/EXAMPLES/small_benchmark_run_to_test_more_complex_Earth/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/EXAMPLES/small_benchmark_run_to_test_very_simple_Earth/DATA/Par_file b/EXAMPLES/small_benchmark_run_to_test_very_simple_Earth/DATA/Par_file index ebd6a5880..eeb42fa5b 100644 --- a/EXAMPLES/small_benchmark_run_to_test_very_simple_Earth/DATA/Par_file +++ b/EXAMPLES/small_benchmark_run_to_test_very_simple_Earth/DATA/Par_file @@ -267,6 +267,7 @@ USE_FORCE_POINT_SOURCE = .false. # half duration is interpreted as a PERIOD just to avoid changing CMTSOLUTION file format # default is .false. which uses Heaviside function USE_MONOCHROMATIC_CMT_SOURCE = .false. +USE_SINSQ_STF = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .false. diff --git a/doc/USER_MANUAL/09_gravity_calculations.tex b/doc/USER_MANUAL/09_gravity_calculations.tex index c57259655..28ccedb71 100644 --- a/doc/USER_MANUAL/09_gravity_calculations.tex +++ b/doc/USER_MANUAL/09_gravity_calculations.tex @@ -1,8 +1,11 @@ -\chapter{Gravity integral calculations for the gravity field of the Earth} +\chapter{Computation of static and transient gravity field} -SPECFEM3D Globe can now compute the gravity field as well as its derivatives (i.e., gravity gradiometry) -generated by any given 3D Earth model at the height +\section{Static gravity field} + +SPECFEM can now compute the static (background) gravity field as well as +its derivatives (i.e., gravity gradiometry) +generated by any given 3D Earth model at the height of an observation satellite, for instance GOCE\newline (see e.g. en.wikipedia.org/wiki/Gravity\_Field\_and\_Steady-State\_Ocean\_Circulation\_Explorer).\newline That feature is still experimental but should work just fine.\newline @@ -16,3 +19,138 @@ \chapter{Gravity integral calculations for the gravity field of the Earth} All main gravity field computations can be found in file \texttt{src/meshfem3d/gravity\_integrals.F90}. Please make sure you compile the code with double-precision, i.e., use flag \texttt{{-}{-}enable-double-precision} for the configuration of the package. + +\section{Coupled elasto-gravitational (transient) simulations} + +Ground deformation modelled in SPECFEM causes motion of +Earth's density field,which produces perturbations in its +gravitational field. This perturbation in the gravitational +field acts as a separate force in the elastodynamics +equations, which introduces itself drives ground deformation. + Hence there is implicit coupling between the deformation + of the ground and the changing gravity field. SPECFEM3D\_GLOBE + typically models gravity using \textit{Cowling's Approximation} + in which only the background gravity field is incorporated, + while the transient perturbation in the gravity field is + ignored. In some cases, however, you may wish to include the + full elastogravitational coupling. This section describes + how to incorporate this coupling into your simulation. + Those using this feature should cite \cite{GhEaTr23}. + + + \subsection{Should I use full gravitational coupling in my simulation?} + + Simulations incorporating full elastogravitational are \textbf{substantially} slower than when using \textit{Cowling's Approximation}. + It is therefore recommended that you only use full gravitational coupling if you really need it. There are two main cases + for which you may be interested in switching this feature on: + + \begin{enumerate} + + \item \textit{ Long-period wavefields} + + The effect of gravity perturbations on the elastic wavefield are negligible at periods below ~ 50-100 seconds; see, for + example, \cite{GhEaTr23} Fig. 13 for a demonstration of the gravitational effect on a wavefield filtered to various + frequencies. Therefore you can ignore these effects unless you are interested in long-period seismology. + + \item \textit{ Deformation-induced gravity modelling} + + Evidently, if you are interested in simulating the gravitational signals associated with ground motion then it is critical + to include self-gravitation! + + \item \textit{ Accurate modelling of an accelerometer} + + As disussed in \cite{DaTr98}, Ch. 4, an accelerometer does not simply measure the ground acceleration. The measured signal + is a combination of the ground acceleration, perturbation in gravitational acceleration, Coriolis effects, and a free-air effect. + Hence, for accurate modelling of acceleration measured at a receiver, one may wish to compute these other terms. These + effects are, for example, very important in the computation of Prompt Elastogravitational Signals (PEGS). Each of these aforementioned + effects measured by an accelerometer can be outputted from SPECFEM3D\_GLOBE. + +\end{enumerate} + + + +\subsection{Requirements \& configuration for the solver} + +Incorporating poisson's equation requires solution to a global matrix of the form $\mathbf{A}\mathbf{x} = \mathbf{b}$. +For this two options are available: an in-built solver, as well as a solver built using the \href{http://petsc.org/release/}{PETSC} toolkit. +Although a built-in solver is implemeted, we recommend using the PETSC solver as it has been benchmarked. As a result, +you will require a local version of PETSC to be installed, which is not shipped with the main SPECFEM3D\_GLOBE package + + +\subsubsection{Installing PETSC} + + PETSC may be installed using a variety of methods depending on your hardware and software. It is recommended that you follow + the \href{https://petsc.org/release/install/}{installation guidelines} on their website, but a few tips are listed below: + + \begin{itemize} + \item PETSC can be built without a fortran compiler for use in c++. Remember to provide a fortran compiler using \texttt{--with-fc} + \item If you are using OpenMPI, you will need to explicitly add the OpenMPI library to your path before configuring. For example using, + \begin{verbatim} + $LD_LIBRARY_PATH=/path/to/OpenMPI/lib/:$LD_LIBRARY_PATH + \end{verbatim} + \item We recommend installing into a separate directory from the source directory by using the \texttt{--prefix} flag + \end{itemize} + + + An example configuration of PETSC for SPECFEM can be found on a \href{here}{https://github.com/williameaton/SpecfemMagic} + on a fork of the \href{SpecfemMagic}{https://github.com/lsawade/SpecfemMagic} scripts originally written by Lucas Sawade. + + + \subsubsection{Compiling SPECFEM3D\_GLOBE with PETSC} + + Compiling SPECFEM3D\_GLOBE with PETSC can be done using the \texttt{--with-petsc} flag. You will also need to include the + path to the PETSC \texttt{include} and \text{lib} directories, as well as add the PETSC library path (PETSC\_LIB) to your + dynamic link library path (LD\_LIBRARY\_PATH). For example + \begin{verbatim} + # define the directory for petsc install + $PETSC_DIR=/path/to/petsc_install_directory + + # petsc library directory + $PETSC_LIB=$PETSC_DIR/lib + + # petsc include library + $PETSC_INC=$PETSC_DIR/include + + # add library to LD library paths + export LD_LIBRARY_PATH=$PETSC_LIB/:$LD_LIBRARY_PATH + + # compile specfem with PETSC: + ./configure CC=$CC CXX=$CXX FC=$FC MPIFC=$MPIFC --with-petsc \ + PETSC_INC=$PETSC_INC PETSC_LIB=$PETSC_LIB + \end{verbatim} + + + \subsection{Running a full-gravity simulation} + A few three changes are required to run with full gravity. First, the \texttt{Par\_file} must be updated as follows: + \begin{enumerate} + \item Both the flags \texttt{GRAVITY} and \texttt{FULL\_GRAVITY} must be set to \texttt{.true.}. + \item Set \texttt{POISSON\_SOLVER} to 0 (in-built solver) or 1 (PETSC solver). PETSC is recommended. + \item Since the current implementation is restricted to global simulations, \texttt{NCHUNKS} must equal 6. + \end{enumerate} +The mesher can then be run as normal. An intermediate workflow step is required between meshing and solving. +Following execution of Meshfem, you must run the program \texttt{xgindex3D}, a serial executable that reorders +some of the global indexing for the new mesh. Simply \texttt{make xgindex3D}, followed by \texttt{./bin/xgindex3D }, +where nprocs is the number of cores/tasks you will use for an MPI run. In some cases you may need to use +\texttt{mpirun -np 1 ./bin/xgindex3D }. Following this, the solver may be executed. + + \subsection{Outputs} + In addition to the normal displacement/velocity/acceleration time-series for a station, a number of other time-series + will be written to \texttt{OUTPUT\_FILES}: + + \begin{itemize} + \item \texttt{net.sta.MXG.sem.ascii} - Gravitational potential + \item \texttt{net.sta.MXZ.C.PGRAV.sem.ascii} - Z direction gravitational acceleration + \end{itemize} + + As discussed above, an idealised vertical accelerometer has four terms: + \begin{equation} + A_v = \mathbf{\hat{\nu}} \cdot \left( \partial_t^2 \mathbf{s} + 2\mathbf{\Omega} \times \partial_t\mathbf{s} + \nabla\phi \right) + \mathbf{s}\cdot \nabla\nabla\Phi \label{eq:idealised_accelerometer} + \end{equation} + for which $\hat{\nu}$ is the direction vector, $\mathbf{s}$ is displacement, ${\Omega}$ is the angular velocity, $\phi$ is the + perturbation in gravitational potential and $\Phi$ is the background gravitational potential. The terms related to Coriolis + and the background potential are then outputted as + + \begin{itemize} + \item \texttt{net.sta.MXZ.C.CORIO.sem.ascii} - Coriolis term + \item \texttt{net.sta.MXZ.C.GRAV.sem.ascii} - Free air term + \end{itemize} \ No newline at end of file diff --git a/doc/USER_MANUAL/bibliography.bib b/doc/USER_MANUAL/bibliography.bib index f2a4a05a5..827c8ae55 100644 --- a/doc/USER_MANUAL/bibliography.bib +++ b/doc/USER_MANUAL/bibliography.bib @@ -8106,6 +8106,17 @@ @article{dehoop2 volume = {64}, year = {1999}} +@article{GhEaTr23, + author = {Gharti, Hom Nath and Eaton, Will and Tromp, Jeroen}, + title = "{Spectral-infinite-element simulations of seismic wave propagation in self-gravitating, rotating 3{D} Earth models}", + journal = {Geophysical Journal International}, + pages = {ggad382}, + year = {2023}, + issn = {0956-540X}, + doi = {10.1093/gji/ggad382}, +} + + @article{HoraCervena2012, author = {Hora, Petr and {\v{C}}erven{\'a}, Olga}, journal = {Applied and Computational Mechanics}, diff --git a/src/shared/broadcast_computed_parameters.f90 b/src/shared/broadcast_computed_parameters.f90 index 045db8970..ce7f29920 100644 --- a/src/shared/broadcast_computed_parameters.f90 +++ b/src/shared/broadcast_computed_parameters.f90 @@ -37,7 +37,7 @@ subroutine broadcast_computed_parameters() integer, parameter :: nparam_i = 51 integer, dimension(nparam_i) :: bcast_integer - integer, parameter :: nparam_l = 80 + integer, parameter :: nparam_l = 81 logical, dimension(nparam_l) :: bcast_logical integer, parameter :: nparam_dp = 42 @@ -113,7 +113,7 @@ subroutine broadcast_computed_parameters() OUTPUT_SEISMOS_3D_ARRAY, & REGIONAL_MESH_CUTOFF,REGIONAL_MESH_ADD_2ND_DOUBLING, & EMC_MODEL,EMC_MODEL_TISO,EMC_MODEL_QMU, & - FULL_GRAVITY, & + FULL_GRAVITY, USE_SINSQ_STF, & HDF5_ENABLED, HDF5_FOR_MOVIES, OUTPUT_SEISMOS_HDF5, & ATTENUATION_3D_BERKELEY /) @@ -367,10 +367,11 @@ subroutine broadcast_computed_parameters() EMC_MODEL_TISO = bcast_logical(74) EMC_MODEL_QMU = bcast_logical(75) FULL_GRAVITY = bcast_logical(76) - HDF5_ENABLED = bcast_logical(77) - HDF5_FOR_MOVIES = bcast_logical(78) - OUTPUT_SEISMOS_HDF5 = bcast_logical(79) - ATTENUATION_3D_BERKELEY = bcast_logical(80) + USE_SINSQ_STF = bcast_logical(77) + HDF5_ENABLED = bcast_logical(78) + HDF5_FOR_MOVIES = bcast_logical(79) + OUTPUT_SEISMOS_HDF5 = bcast_logical(80) + ATTENUATION_3D_BERKELEY = bcast_logical(81) ! double precisions DT = bcast_double_precision(1) diff --git a/src/shared/read_parameter_file.F90 b/src/shared/read_parameter_file.F90 index 7adaf30b0..efab5452a 100644 --- a/src/shared/read_parameter_file.F90 +++ b/src/shared/read_parameter_file.F90 @@ -216,6 +216,11 @@ subroutine read_parameter_file() call read_value_logical(USE_FORCE_POINT_SOURCE, 'USE_FORCE_POINT_SOURCE', ier) if (ier /= 0) stop 'an error occurred while reading the parameter file: USE_FORCE_POINT_SOURCE' + ! sin squared STF - used for PEGS + call read_value_logical(USE_SINSQ_STF, 'USE_SINSQ_STF', ier) + if (ier /= 0) stop 'an error occurred while reading the parameter file: USE_SINSQ_STF' + + ! monochromatic double couple source call read_value_logical(USE_MONOCHROMATIC_CMT_SOURCE, 'USE_MONOCHROMATIC_CMT_SOURCE', ier) if (ier /= 0) stop 'an error occurred while reading the parameter file: USE_MONOCHROMATIC_CMT_SOURCE' diff --git a/src/shared/shared_par.f90 b/src/shared/shared_par.f90 index 108f64074..cbb62f510 100644 --- a/src/shared/shared_par.f90 +++ b/src/shared/shared_par.f90 @@ -76,6 +76,7 @@ module shared_input_parameters ! sources logical :: USE_FORCE_POINT_SOURCE + logical :: USE_SINSQ_STF = .false. logical :: USE_MONOCHROMATIC_CMT_SOURCE,PRINT_SOURCE_TIME_FUNCTION ! checkpointing/restart diff --git a/src/specfem3D/SIEM_compute_seismograms.F90 b/src/specfem3D/SIEM_compute_seismograms.F90 index 48172ac5b..606080563 100644 --- a/src/specfem3D/SIEM_compute_seismograms.F90 +++ b/src/specfem3D/SIEM_compute_seismograms.F90 @@ -86,7 +86,7 @@ subroutine SIEM_compute_seismograms_phi(phi_crust_mantle,seismograms) use specfem_par_crustmantle, only: ibool_crust_mantle - use specfem_par_full_gravity, only: scale_pgrav + use specfem_par_full_gravity, only: scale_phi implicit none @@ -121,7 +121,7 @@ subroutine SIEM_compute_seismograms_phi(phi_crust_mantle,seismograms) ! store North, East and Vertical components ! distinguish between single and double precision for reals - seismograms(:,irec_local,seismo_current) = real(scale_pgrav * phir,kind=CUSTOM_REAL) + seismograms(:,irec_local,seismo_current) = real(scale_phi * phir,kind=CUSTOM_REAL) enddo end subroutine SIEM_compute_seismograms_phi @@ -147,7 +147,7 @@ subroutine SIEM_compute_seismograms_pgrav(var_crust_mantle, seismograms, seismog etax_crust_mantle, etay_crust_mantle, etaz_crust_mantle, & gammax_crust_mantle, gammay_crust_mantle, gammaz_crust_mantle - use specfem_par_full_gravity, only: scale_pgrav + use specfem_par_full_gravity, only: scale_pgrav, scale_hgrav implicit none @@ -285,7 +285,7 @@ subroutine SIEM_compute_seismograms_pgrav(var_crust_mantle, seismograms, seismog eps_loc_rot(:,:) = real(matmul(nu_rec(:,:,irec), Hgrav_rec(:,:)),kind=CUSTOM_REAL) eps_loc_new(:,:) = real(matmul(eps_loc_rot(:,:), transpose(nu_rec(:,:,irec))),kind=CUSTOM_REAL) - seismograms_Hgrav(:,:,irec_local,seismo_current) = real(scale_pgrav * eps_loc_new(:,:),kind=CUSTOM_REAL) + seismograms_Hgrav(:,:,irec_local,seismo_current) = real(scale_hgrav * eps_loc_new(:,:),kind=CUSTOM_REAL) ! original: misses P' projection for rotation !do m = 1,3 diff --git a/src/specfem3D/SIEM_prepare_iteration.F90 b/src/specfem3D/SIEM_prepare_iteration.F90 index bfd103c09..2c4808f0f 100644 --- a/src/specfem3D/SIEM_prepare_iteration.F90 +++ b/src/specfem3D/SIEM_prepare_iteration.F90 @@ -451,6 +451,7 @@ end subroutine write_ensight_perelementAS call compute_poisson_rhoload3() if (POISSON_SOLVER == ISOLVER_BUILTIN) then + if (myrank == 0) print *,' Using built-in solver...' ! built-in solver if (CG_SCALING) then gravload1(:) = ndscale1(:) * gravload1(:) @@ -465,6 +466,7 @@ end subroutine write_ensight_perelementAS endif else ! petsc solver + if (myrank == 0) print *,' Using PETSC solver...' call petsc_set_vector1(gravload1) if (myrank == 0) print *,'PETSc solver: L1: Set RHS with PETSC done' @@ -763,6 +765,7 @@ end subroutine write_ensight_perelementAS write(IMAIN,*) ' Elapsed time for background gravity solve: ',sngl(tCPU),'(s)' write(IMAIN,*) call flush_IMAIN() + print *, '--- COMPLETED BACKGROUND GRAVITY ---' endif end subroutine SIEM_compute_background_gravity diff --git a/src/specfem3D/SIEM_prepare_solver.F90 b/src/specfem3D/SIEM_prepare_solver.F90 index 25344f4f9..e2b2e0a3f 100644 --- a/src/specfem3D/SIEM_prepare_solver.F90 +++ b/src/specfem3D/SIEM_prepare_solver.F90 @@ -1635,8 +1635,15 @@ subroutine SIEM_prepare_seismos() if (.not. FULL_GRAVITY) return ! additional scaling factors for gravity seismograms - scale_accel = scale_veloc * scale_t_inv ! [m / s^2] - scale_pgrav = scale_displ**2 * scale_t_inv**2 ! [m^2 / s^2] ONE NEED TO BE CHECKED!!! + scale_accel = scale_veloc * scale_t_inv ! [m / s^2] + scale_phi = scale_displ**2 * scale_t_inv**2 ! [m^2 / s^2] + scale_pgrav = scale_displ * scale_t_inv**2 ! [m / s^2] + + ! strain is defined as second time integral + ! of grad grad phi here we just output + ! grad grad phi so scale by units of 1/s^2 + scale_hgrav = scale_t_inv**2 ! [1 / s^2] + ! allocate seismogram array if (SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3) then diff --git a/src/specfem3D/SIEM_solver_petsc.F90 b/src/specfem3D/SIEM_solver_petsc.F90 index cc60fc891..d0aaa0849 100644 --- a/src/specfem3D/SIEM_solver_petsc.F90 +++ b/src/specfem3D/SIEM_solver_petsc.F90 @@ -80,6 +80,24 @@ ! call VecSetValues(bvec,neq,l2gdof(1:),y,ADD_VALUES,ierr) ! well, until a better solution is found. +! 5/12/25 - WE: A further issue with this is the variables ix and y +! for single values e.g. originally +! call VecSetValues(interface_gvec1,1,igdof,rval,INSERT_VALUES,ierr) +! needs to be written as +! call VecSetValues(interface_gvec1,1,[igdof],[rval],INSERT_VALUES,ierr) +! also for MatSetValues(Mat mat, PetscInt m, const PetscInt idxm[], +! PetscInt n, const PetscInt idxn[], +! const PetscScalar v[], InsertMode addv) +! the variables idxm, idxn, and v. +! Other examples in: +! MatMPIAIJSetPreallocation --> requires PETSC_NULL_INTEGER_ARRAY +! instead of PETSC_NULL_INTEGER + +! For MatSetValues - flattening +! MatSetValues documentation suggests in f90 varray needs to be flattened +! note that F90 is column major but C is row major. +! It appears petsc updated VecGetArrayF90 back to VecGetArray in +! 3.23 (https://petsc.org/release/changes/323/) module siem_solver_petsc @@ -117,7 +135,6 @@ module siem_solver_petsc !use petscksp, only: tKSP !use petscpc, only: tPC !use petscis, only: tIS - #endif use constants, only: myrank,IMAIN,CUSTOM_REAL @@ -178,17 +195,48 @@ module siem_solver_petsc ! KSPSetTolerances(KSP ksp, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt maxits) ! ! Level-1 KSP solver - PetscInt, parameter :: KSP_MAXITER1 = 3000 + PetscInt, parameter :: KSP_MAXITER1 = 3000 PetscReal, parameter :: KSP_RTOL1 = 1.0e-7 PetscReal, parameter :: KSP_ATOL1 = 1.0e-30 PetscReal, parameter :: KSP_DTOL1 = 1.0e30 ! Level-2 KSP solver - PetscInt, parameter :: KSP_MAXITER = 3000 + PetscInt, parameter :: KSP_MAXITER = 3000 PetscReal, parameter :: KSP_RTOL = 1.0e-7 PetscReal, parameter :: KSP_ATOL = 1.0e-30 PetscReal, parameter :: KSP_DTOL = 1.0e30 + ! Centralize PETSc version check + ! Force display of version info (will cause compilation to stop here) + ! PETSC_NULL_SCALAR_ARRAY instead of PETSC_NULL_SCALAR + ! introduced in 3.22 +#if PETSC_VERSION_GE(3,22,0) +#define PETSC_ARRAY_NULL_API 1 +#else +#define PETSC_ARRAY_NULL_API 0 +#endif + +! Versions between 3.19 and 3.23 used an F90 +! suffix in things like VecGetArray +#if PETSC_VERSION_GE(3,19,0) +#if PETSC_VERSION_LT(3,23,0) +#define API_F90_SUFFIX 1 +#else +#define API_F90_SUFFIX 0 +#endif +#else +#define API_F90_SUFFIX 0 +#endif + +! In version 3.23 +! MatInfo changed from arrays to derived data types +#if PETSC_VERSION_GE(3,22,0) +#define PETSC_INFO_DERIVED_TYPE 1 +#else +#define PETSC_INFO_DERIVED_TYPE 0 +#endif + + ! solver type PetscInt, parameter :: COMMAND = 0, CG = 1, SUPERLU = 2, MUMPS = 3 @@ -234,7 +282,7 @@ module siem_solver_petsc type(tPC) :: pc PetscErrorCode :: ierr PetscInt :: nzeros_max,nzeros_min - PetscInt :: ig0,ig1 + PetscInt :: ig0,ig1, ONE #endif ! public function @@ -264,7 +312,7 @@ subroutine petsc_initialize1() use specfem_par, only: ADD_TRINF,SIMULATION_TYPE use specfem_par_full_gravity, only: ggdof_ic1,ggdof_oc1,ggdof_cm1,ggdof_trinf1,ggdof_inf1 - + implicit none type(tVec) :: nzeror_gvec1,nzeror_dvec1,nzeror_ovec1,iproc_gvec1, & interface_gvec1,ninterface_dvec1,ninterface_ovec1,nself_gvec1 @@ -293,7 +341,6 @@ subroutine petsc_initialize1() PetscInt,allocatable :: ninterface_darray1(:),ninterface_oarray1(:) PetscScalar,allocatable :: rg_interface(:),rnself_lgarray1(:) PetscScalar,pointer :: rninterface_darray1(:),rninterface_oarray1(:) - ! memory info PetscLogDouble :: bytes @@ -479,7 +526,7 @@ subroutine petsc_initialize1() do i_bool = 1,nibool do i_ndof = 1,NNDOF igdof = ggdof_ic1(i_ndof,ibool_interface(i_bool))-1 - if (igdof >= 0) call VecSetValues(interface_gvec1,1,igdof,rval,INSERT_VALUES,ierr) + if (igdof >= 0) call VecSetValues(interface_gvec1,1,[igdof],[rval],INSERT_VALUES,ierr) enddo enddo deallocate(ibool_interface) @@ -493,7 +540,7 @@ subroutine petsc_initialize1() do i_bool = 1,nibool do i_ndof = 1,NNDOF igdof = ggdof_oc1(i_ndof,ibool_interface(i_bool))-1 - if (igdof >= 0) call VecSetValues(interface_gvec1,1,igdof,rval,INSERT_VALUES,ierr) + if (igdof >= 0) call VecSetValues(interface_gvec1,1,[igdof],[rval],INSERT_VALUES,ierr) enddo enddo deallocate(ibool_interface) @@ -507,7 +554,7 @@ subroutine petsc_initialize1() do i_bool = 1,nibool do i_ndof = 1,NNDOF igdof = ggdof_cm1(i_ndof,ibool_interface(i_bool))-1 - if (igdof >= 0) call VecSetValues(interface_gvec1,1,igdof,rval,INSERT_VALUES,ierr) + if (igdof >= 0) call VecSetValues(interface_gvec1,1,[igdof],[rval],INSERT_VALUES,ierr) enddo enddo deallocate(ibool_interface) @@ -522,7 +569,7 @@ subroutine petsc_initialize1() do i_bool = 1,nibool do i_ndof = 1,NNDOF igdof = ggdof_trinf1(i_ndof,ibool_interface(i_bool))-1 - if (igdof >= 0) call VecSetValues(interface_gvec1,1,igdof,rval,INSERT_VALUES,ierr) + if (igdof >= 0) call VecSetValues(interface_gvec1,1,[igdof],[rval],INSERT_VALUES,ierr) enddo enddo deallocate(ibool_interface) @@ -537,7 +584,7 @@ subroutine petsc_initialize1() do i_bool = 1,nibool do i_ndof = 1,NNDOF igdof = ggdof_inf1(i_ndof,ibool_interface(i_bool))-1 - if (igdof >= 0) call VecSetValues(interface_gvec1,1,igdof,rval,INSERT_VALUES,ierr) + if (igdof >= 0) call VecSetValues(interface_gvec1,1,[igdof],[rval],INSERT_VALUES,ierr) enddo enddo deallocate(ibool_interface) @@ -559,7 +606,7 @@ subroutine petsc_initialize1() rval = 1.0 do i = 1,neq1 if (isg_interface(i) == 1) then - call VecSetValues(nself_gvec1,1,l2gdof1(i),rval,ADD_VALUES,ierr); + call VecSetValues(nself_gvec1,1,[l2gdof1(i)],[rval],ADD_VALUES,ierr); endif enddo call VecAssemblyBegin(nself_gvec1,ierr); CHECK_PETSC_ERROR(ierr) @@ -568,14 +615,24 @@ subroutine petsc_initialize1() allocate(rnself_lgarray1(neq1)) call scatter_globalvec1(nself_gvec1, rnself_lgarray1) + +#if API_F90_SUFFIX call VecGetArrayF90(nself_gvec1,rnself_array1,ierr) +#else + call VecGetArray(nself_gvec1,rnself_array1,ierr) +#endif allocate(nself_array1(n)) nself_array1 = int(rnself_array1(1:n)) where(nself_array1 > 0) nself_array1 = nself_array1-1 ! subtract self +#if API_F90_SUFFIX call VecRestoreArrayF90(nself_gvec1,rnself_array1,ierr) +#else + call VecRestoreArray(nself_gvec1,rnself_array1,ierr) +#endif + call VecDestroy(nself_gvec1,ierr) if (myrank == 0) print *,'PETSc solver: maximum value of nself:',maxval(nself_array1) @@ -619,11 +676,11 @@ subroutine petsc_initialize1() ! set values computed so far rnd = real(nd) rnoffd = real(noffd) - call VecSetValues(nzeror_dvec1,1,igr0,rnd,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) - call VecSetValues(nzeror_ovec1,1,igr0,rnoffd,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call VecSetValues(nzeror_dvec1,1,[igr0],[rnd],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call VecSetValues(nzeror_ovec1,1,[igr0],[rnoffd],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) - call VecSetValues(ninterface_dvec1,1,igr0,rnid,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) - call VecSetValues(ninterface_ovec1,1,igr0,rnioffd,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call VecSetValues(ninterface_dvec1,1,[igr0],[rnid],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call VecSetValues(ninterface_ovec1,1,[igr0],[rnioffd],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) ! reset nd = 0; noffd = 0 @@ -667,11 +724,11 @@ subroutine petsc_initialize1() ! for last rnd = real(nd) rnoffd = real(noffd) - call VecSetValues(nzeror_dvec1,1,igr0,rnd,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) - call VecSetValues(nzeror_ovec1,1,igr0,rnoffd,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call VecSetValues(nzeror_dvec1,1,[igr0],[rnd],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call VecSetValues(nzeror_ovec1,1,[igr0],[rnoffd],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) - call VecSetValues(ninterface_dvec1,1,igr0,rnid,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) - call VecSetValues(ninterface_ovec1,1,igr0,rnioffd,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call VecSetValues(ninterface_dvec1,1,[igr0],[rnid],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call VecSetValues(ninterface_ovec1,1,[igr0],[rnioffd],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) endif enddo deallocate(krow_sparse1,kcol_sparse1) @@ -690,7 +747,12 @@ subroutine petsc_initialize1() ! apply correction for repeatition due to interfaces ! diagonal matrix call VecGetLocalSize(nzeror_dvec1,n,ierr) +#if API_F90_SUFFIX call VecGetArrayF90(nzeror_dvec1,nzeror_darray1,ierr) +#else + call VecGetArray(nzeror_dvec1,nzeror_darray1,ierr) +#endif + allocate(nnzero_diag1(n)) nnzero_diag1(:) = int(nzeror_darray1(1:n)) @@ -701,42 +763,77 @@ subroutine petsc_initialize1() minval(nnzero_diag1),maxval(nnzero_diag1) call synchronize_all() +#if API_F90_SUFFIX call VecRestoreArrayF90(nzeror_dvec1,nzeror_darray1,ierr) +#else + call VecRestoreArray(nzeror_dvec1,nzeror_darray1,ierr) +#endif + call VecDestroy(nzeror_dvec1,ierr) ! off-diagonal matrix +#if API_F90_SUFFIX call VecGetArrayF90(nzeror_ovec1,nzeror_oarray1,ierr) +#else + call VecGetArray(nzeror_ovec1,nzeror_oarray1,ierr) +#endif + allocate(nnzero_offdiag1(n)) nnzero_offdiag1(:) = int(nzeror_oarray1(1:n)) +#if API_F90_SUFFIX call VecRestoreArrayF90(nzeror_ovec1,nzeror_oarray1,ierr) +#else + call VecRestoreArray(nzeror_ovec1,nzeror_oarray1,ierr) +#endif + call VecDestroy(nzeror_ovec1,ierr) ! correction ! I do not know why but there are some DOFs where the correction exceeds by 4 or ! 8 therefore to be safe we need to subtract this from all +#if API_F90_SUFFIX call VecGetArrayF90(ninterface_dvec1,rninterface_darray1,ierr) +#else + call VecGetArray(ninterface_dvec1,rninterface_darray1,ierr) +#endif + !where(rninterface_darray1>0.0 .and. rninterface_darray1 < 1.0)rninterface_darray1=1.0 allocate(ninterface_darray1(n)) ninterface_darray1 = int(rninterface_darray1(1:n)) +#if API_F90_SUFFIX call VecRestoreArrayF90(ninterface_dvec1,rninterface_darray1,ierr) +#else + call VecRestoreArray(ninterface_dvec1,rninterface_darray1,ierr) +#endif + call VecDestroy(ninterface_dvec1,ierr) where(ninterface_darray1 > 0) ninterface_darray1 = ninterface_darray1 - 4 where(ninterface_darray1 < 0) ninterface_darray1 = 0 +#if API_F90_SUFFIX call VecGetArrayF90(ninterface_ovec1,rninterface_oarray1,ierr) +#else + call VecGetArray(ninterface_ovec1,rninterface_oarray1,ierr) +#endif + !where(rninterface_oarray1>0.0 .and. rninterface_oarray1 < 1.0)rninterface_oarray1=1.0 allocate(ninterface_oarray1(n)) ninterface_oarray1 = int(rninterface_oarray1(1:n)) - call VecRestoreArrayF90(ninterface_ovec1,rninterface_oarray1,ierr) +#if API_F90_SUFFIX + call VecGetArrayF90(ninterface_ovec1,rninterface_oarray1,ierr) +#else + call VecGetArray(ninterface_ovec1,rninterface_oarray1,ierr) +#endif + call VecDestroy(ninterface_ovec1,ierr) where(ninterface_oarray1 > 0) ninterface_oarray1 = ninterface_oarray1 - 8 @@ -753,7 +850,7 @@ subroutine petsc_initialize1() rval = 1.0 do i = 1,nsparse1 igdof = kgrow_sparse1(i)-1 ! Fortran index - call VecSetValues(nzeror_gvec1,1,igdof,rval,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call VecSetValues(nzeror_gvec1,1,[igdof],[rval],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) enddo call VecAssemblyBegin(nzeror_gvec1,ierr); CHECK_PETSC_ERROR(ierr) call VecAssemblyEnd(nzeror_gvec1,ierr); CHECK_PETSC_ERROR(ierr) @@ -768,12 +865,21 @@ subroutine petsc_initialize1() deallocate(kgrow_sparse1,kgcol_sparse1) ! non-zero array for diagonal/off-diagonal matrix? +#if API_F90_SUFFIX call VecGetArrayF90(nzeror_gvec1,nzeror_array1,ierr); CHECK_PETSC_ERROR(ierr) +#else + call VecGetArray(nzeror_gvec1,nzeror_array1,ierr); CHECK_PETSC_ERROR(ierr) +#endif allocate(inzeror_array1(n)) inzeror_array1(:) = int(nzeror_array1(1:n)) +#if API_F90_SUFFIX call VecRestoreArrayF90(nzeror_gvec1,nzeror_array1,ierr); CHECK_PETSC_ERROR(ierr) +#else + call VecRestoreArray(nzeror_gvec1,nzeror_array1,ierr); CHECK_PETSC_ERROR(ierr) +#endif + call VecDestroy(nzeror_gvec1,ierr); CHECK_PETSC_ERROR(ierr) inzeros_max = maxvec(inzeror_array1) @@ -808,8 +914,19 @@ subroutine petsc_initialize1() ! ! this seems to lead to a much faster petsc_set_matrix1() routine without the re-allocations. ! however, the diagonal and in particular the off-diagonal estimate with nzeros_max might be still off. - call MatMPIAIJSetPreallocation(Amat1,nzeros_max,PETSC_NULL_INTEGER, & - nzeros_max,PETSC_NULL_INTEGER,ierr); CHECK_PETSC_ERROR(ierr) + + +#if PETSC_ARRAY_NULL_API + ! Later version uses integer array + call MatMPIAIJSetPreallocation(Amat1,nzeros_max, PETSC_NULL_INTEGER_ARRAY, & + nzeros_max, PETSC_NULL_INTEGER_ARRAY, ierr); CHECK_PETSC_ERROR(ierr) +#else + call MatMPIAIJSetPreallocation(Amat1,nzeros_max, PETSC_NULL_INTEGER, & + nzeros_max, PETSC_NULL_INTEGER, ierr); CHECK_PETSC_ERROR(ierr) +#endif + + + call MatSetFromOptions(Amat1,ierr); CHECK_PETSC_ERROR(ierr) call MatGetOwnershipRange(Amat1,istart,iend,ierr); CHECK_PETSC_ERROR(ierr) @@ -1077,8 +1194,13 @@ subroutine petsc_set_matrix1() !debugging logical, parameter :: DEBUG_FILE_OUTPUT = .false. - integer :: ncols - integer,dimension(:),allocatable :: cols + + ! For the MatGetRow call + !integer :: ncols + !integer,dimension(:),allocatable :: cols + PetscInt :: ncols + PetscInt, pointer :: cols(:) + PetscScalar, pointer :: debug_vals(:) character(len=10) :: char_myrank character(len=60) :: outf_name @@ -1088,9 +1210,16 @@ subroutine petsc_set_matrix1() ! const PetscInt idxn[], const PetscScalar v[], InsertMode addv) PetscScalar :: v - ! matrix info + ! matrix info +#if PETSC_INFO_DERIVED_TYPE +MatInfo :: info +#else double precision :: info(MAT_INFO_SIZE) - double precision :: mallocsval +#endif +double precision :: mallocsval + + + PetscLogDouble :: bytes ! timing @@ -1150,7 +1279,7 @@ subroutine petsc_set_matrix1() ! storekmat_inner_core1(i,j,i_elmt),ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) ! petsc types v = storekmat_inner_core1(i,j,i_elmt) - call MatSetValues(Amat1,1,ggdof_elmt(i),1,ggdof_elmt(j),v,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call MatSetValues(Amat1,1,[ggdof_elmt(i)],1,[ggdof_elmt(j)],[v],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) endif enddo enddo @@ -1179,7 +1308,7 @@ subroutine petsc_set_matrix1() ! storekmat_outer_core1(i,j,i_elmt),ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) ! petsc types v = storekmat_outer_core1(i,j,i_elmt) - call MatSetValues(Amat1,1,ggdof_elmt(i),1,ggdof_elmt(j),v,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call MatSetValues(Amat1,1,[ggdof_elmt(i)],1,[ggdof_elmt(j)],[v],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) endif enddo enddo @@ -1208,7 +1337,7 @@ subroutine petsc_set_matrix1() ! storekmat_crust_mantle1(i,j,i_elmt),ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) ! petsc types v = storekmat_crust_mantle1(i,j,i_elmt) - call MatSetValues(Amat1,1,ggdof_elmt(i),1,ggdof_elmt(j),v,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call MatSetValues(Amat1,1,[ggdof_elmt(i)],1,[ggdof_elmt(j)],[v],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) endif enddo enddo @@ -1243,7 +1372,7 @@ subroutine petsc_set_matrix1() ! storekmat_trinfinite1(i,j,i_elmt),ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) ! petsc types v = storekmat_trinfinite1(i,j,i_elmt) - call MatSetValues(Amat1,1,ggdof_elmt(i),1,ggdof_elmt(j),v,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call MatSetValues(Amat1,1,[ggdof_elmt(i)],1,[ggdof_elmt(j)],[v],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) endif enddo enddo @@ -1276,7 +1405,7 @@ subroutine petsc_set_matrix1() ! storekmat_infinite1(i,j,i_elmt),ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) ! petsc types v = storekmat_infinite1(i,j,i_elmt) - call MatSetValues(Amat1,1,ggdof_elmt(i),1,ggdof_elmt(j),v,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call MatSetValues(Amat1,1,[ggdof_elmt(i)],1,[ggdof_elmt(j)],[v],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) endif enddo enddo @@ -1297,26 +1426,29 @@ subroutine petsc_set_matrix1() ! symmetric call MatSetOption(Amat1,MAT_SYMMETRIC,PETSC_TRUE,ierr); CHECK_PETSC_ERROR(ierr) - ! debugging output + !debugging output + !unsure if this updated debugging call works but it compiles... + !updates require pointers if (DEBUG_FILE_OUTPUT) then - allocate(cols(nzeros_max)) - cols(:) = 0 + ! allocate(cols(nzeros_max)) + ! cols(:) = 0 + write(char_myrank,'(i4)') myrank outf_name='tmp_nonzeros'//trim(adjustl(char_myrank)) open(1,file=outf_name,action='write',status='replace') call MatGetOwnershipRange(Amat1,istart,iend,ierr); CHECK_PETSC_ERROR(ierr) do i = istart,iend-1 - cols(:) = -1 + !cols(:) = -1 ! gets row i - call MatGetRow(Amat1,i,ncols,cols,PETSC_NULL_SCALAR,ierr); CHECK_PETSC_ERROR(ierr) + call MatGetRow(Amat1, i, ncols, cols, debug_vals,ierr); CHECK_PETSC_ERROR(ierr) ndiag = count(cols >= ig0 .and. cols <= ig1) noffdiag = ncols-ndiag write(1,*) ndiag,noffdiag,ncols ! free temporary space of MatGetRow() - call MatRestoreRow(Amat1,i,ncols,PETSC_NULL_INTEGER,PETSC_NULL_SCALAR,ierr); CHECK_PETSC_ERROR(ierr) + call MatRestoreRow(Amat1,i,ncols,cols,debug_vals,ierr); CHECK_PETSC_ERROR(ierr) enddo close(1) - deallocate(cols) + !deallocate(cols) endif ! synchronize all processes @@ -1325,7 +1457,11 @@ subroutine petsc_set_matrix1() ! matrix info call MatGetInfo(Amat1, MAT_GLOBAL_MAX, info, ierr); CHECK_PETSC_ERROR(ierr) - mallocsval = info(MAT_INFO_MALLOCS) ! number of mallocs during MatSetValues() +#if PETSC_INFO_DERIVED_TYPE +mallocsval = info%mallocs +#else +mallocsval = info(MAT_INFO_MALLOCS) ! number of mallocs during MatSetValues() +#endif !memval = info(MAT_INFO_MEMORY) ! memory allocated - not provided !nonzeros_allocated = info(MAT_INFO_NZ_ALLOCATED) ! nonzero entries allocated @@ -1474,7 +1610,8 @@ subroutine petsc_solve1(sdata1,niter) #ifdef USE_PETSC ! local parameters - PetscInt :: iter,ireason + PetscInt :: iter + KSPConvergedReason:: ireason ! petsc type array PetscScalar :: y(size(sdata1)) @@ -1612,10 +1749,20 @@ subroutine scatter_globalvec1(global_vec,larray) call VecScatterEnd(vscat1,global_vec,local_vec1,INSERT_VALUES,SCATTER_FORWARD,ierr); CHECK_PETSC_ERROR(ierr) call VecGetSize(local_vec1,n,ierr) +#if API_F90_SUFFIX call VecGetArrayF90(local_vec1,array_data,ierr); CHECK_PETSC_ERROR(ierr) +#else + call VecGetArray(local_vec1,array_data,ierr); CHECK_PETSC_ERROR(ierr) +#endif + larray(1:n) = array_data(1:n) + +#if API_F90_SUFFIX call VecRestoreArrayF90(local_vec1,array_data,ierr); CHECK_PETSC_ERROR(ierr) +#else + call VecRestoreArray(local_vec1,array_data,ierr); CHECK_PETSC_ERROR(ierr) +#endif end subroutine scatter_globalvec1 @@ -1860,7 +2007,7 @@ subroutine petsc_initialize() call synchronize_all() ! preallocation - call MatMPIAIJSetPreallocation(Amat,nzeros_max,nzeros,nzeros_max,20*nzeros,ierr); CHECK_PETSC_ERROR(ierr) + call MatMPIAIJSetPreallocation(Amat,nzeros_max,[nzeros],nzeros_max,[20*nzeros],ierr); CHECK_PETSC_ERROR(ierr) call MatSetFromOptions(Amat,ierr); CHECK_PETSC_ERROR(ierr) call MatGetOwnershipRange(Amat,istart,iend,ierr); CHECK_PETSC_ERROR(ierr) @@ -1981,18 +2128,23 @@ subroutine petsc_set_matrix() ggdof_inf,storekmat_infinite,inode_elmt_inf implicit none - integer :: i,i_elmt,j,ncount + integer :: i,i_elmt,j,ncount,iflat,jflat integer :: ggdof_elmt(NEDOF),idof(NEDOF),igdof(NEDOF) ! types required by MatSetValues: ! MatSetValues(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, \ ! const PetscInt idxn[], const PetscScalar v[], InsertMode addv) PetscScalar :: v - PetscScalar,dimension(:,:),allocatable :: varr + PetscScalar, dimension(:),allocatable :: varr ! matrix info +#if PETSC_INFO_DERIVED_TYPE +MatInfo :: info +#else double precision :: info(MAT_INFO_SIZE) - double precision :: mallocsval +#endif +double precision :: mallocsval + PetscLogDouble :: bytes ! timing @@ -2036,7 +2188,7 @@ subroutine petsc_set_matrix() ! storekmat_inner_core(i,j,i_elmt),ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) ! petsc types v = storekmat_inner_core(i,j,i_elmt) - call MatSetValues(Amat,1,ggdof_elmt(i),1,ggdof_elmt(j),v,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call MatSetValues(Amat,1,[ggdof_elmt(i)],1,[ggdof_elmt(j)],[v],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) endif enddo enddo @@ -2062,7 +2214,7 @@ subroutine petsc_set_matrix() storekmat_outer_core(i,j,i_elmt) /= 0.0_CUSTOM_REAL) then ! petsc types v = storekmat_outer_core(i,j,i_elmt) - call MatSetValues(Amat,1,ggdof_elmt(i),1,ggdof_elmt(j),v,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) + call MatSetValues(Amat,1,[ggdof_elmt(i)],1,[ggdof_elmt(j)],[v],ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) endif enddo enddo @@ -2092,8 +2244,17 @@ subroutine petsc_set_matrix() !call MatSetValues(Amat,ncount,igdof(1:ncount),ncount,igdof(1:ncount), & ! storekmat_crust_mantle(idof(1:ncount),idof(1:ncount),i_elmt),ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) ! petsc types - allocate(varr(ncount,ncount)) - varr(:,:) = storekmat_crust_mantle(idof(1:ncount),idof(1:ncount),i_elmt) + allocate(varr(ncount*ncount)) + + ! WE 5/12/25 + ! See notes at top on flattening + !varr(:,:) = storekmat_crust_mantle(idof(1:ncount),idof(1:ncount),i_elmt) + do iflat = 1, ncount + do jflat = 1, ncount + varr((iflat-1)*ncount + jflat) = storekmat_crust_mantle(idof(iflat), idof(jflat), i_elmt) + enddo !jflat + enddo !iflat + call MatSetValues(Amat,ncount,igdof(1:ncount),ncount,igdof(1:ncount),varr,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) deallocate(varr) enddo @@ -2128,8 +2289,14 @@ subroutine petsc_set_matrix() !call MatSetValues(Amat,ncount,igdof(1:ncount),ncount,igdof(1:ncount), & ! storekmat_trinfinite(idof(1:ncount),idof(1:ncount),i_elmt),ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) ! petsc types - allocate(varr(ncount,ncount)) - varr(:,:) = storekmat_trinfinite(idof(1:ncount),idof(1:ncount),i_elmt) + allocate(varr(ncount*ncount)) + ! WE 5/12/25 + ! See notes at top on flattening + do iflat = 1, ncount + do jflat = 1, ncount + varr((iflat-1)*ncount + jflat) = storekmat_trinfinite(idof(iflat), idof(jflat), i_elmt) + enddo !jflat + enddo !iflat call MatSetValues(Amat,ncount,igdof(1:ncount),ncount,igdof(1:ncount),varr,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) deallocate(varr) enddo @@ -2163,8 +2330,12 @@ subroutine petsc_set_matrix() !call MatSetValues(Amat,ncount,igdof(1:ncount),ncount,igdof(1:ncount), & ! storekmat_infinite(idof(1:ncount),idof(1:ncount),i_elmt),ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) ! petsc types - allocate(varr(ncount,ncount)) - varr(:,:) = storekmat_infinite(idof(1:ncount),idof(1:ncount),i_elmt) + allocate(varr(ncount*ncount)) + do iflat = 1, ncount + do jflat = 1, ncount + varr((iflat-1)*ncount + jflat) = storekmat_infinite(idof(iflat), idof(jflat), i_elmt) + enddo !jflat + enddo !iflat call MatSetValues(Amat,ncount,igdof(1:ncount),ncount,igdof(1:ncount),varr,ADD_VALUES,ierr); CHECK_PETSC_ERROR(ierr) deallocate(varr) enddo @@ -2183,7 +2354,12 @@ subroutine petsc_set_matrix() ! matrix info call MatGetInfo(Amat, MAT_GLOBAL_MAX, info, ierr); CHECK_PETSC_ERROR(ierr) - mallocsval = info(MAT_INFO_MALLOCS) ! number of mallocs during MatSetValues() + +#if PETSC_INFO_DERIVED_TYPE +mallocsval = info%mallocs +#else +mallocsval = info(MAT_INFO_MALLOCS) ! number of mallocs during MatSetValues() +#endif ! memory usage call PetscMemoryGetCurrentUsage(bytes, ierr); CHECK_PETSC_ERROR(ierr) @@ -2282,7 +2458,8 @@ subroutine petsc_solve(sdata,niter) #ifdef USE_PETSC ! local parameters - PetscInt :: iter,ireason + PetscInt :: iter + KSPConvergedReason:: ireason ! petsc type array PetscScalar :: y(size(sdata)) @@ -2295,7 +2472,7 @@ subroutine petsc_solve(sdata,niter) ! Check solution and clean up !------------------------------------------------------------------------------- - call KSPGetConvergedReason(ksp,ireason,ierr); CHECK_PETSC_ERROR(ierr) + call KSPGetConvergedReason(ksp, ireason,ierr); CHECK_PETSC_ERROR(ierr) call KSPGetIterationNumber(ksp,iter,ierr); CHECK_PETSC_ERROR(ierr) !debug diff --git a/src/specfem3D/comp_source_time_function.f90 b/src/specfem3D/comp_source_time_function.f90 index 6f5c11e44..b6f4e373e 100644 --- a/src/specfem3D/comp_source_time_function.f90 +++ b/src/specfem3D/comp_source_time_function.f90 @@ -32,7 +32,7 @@ double precision function comp_source_time_function(t,hdur,it_index) ! for berkeley source time function use shared_parameters, only: STF_IS_UCB_HEAVISIDE use ucb_heaviside, only: comp_source_time_function_ucb_stf - + use specfem_par, only: USE_SINSQ_STF implicit none double precision,intent(in) :: t,hdur @@ -40,6 +40,7 @@ double precision function comp_source_time_function(t,hdur,it_index) ! local parameters double precision, external :: comp_source_time_function_heavi + double precision, external :: comp_source_time_function_ssq double precision, external :: comp_source_time_function_ext if (STF_IS_UCB_HEAVISIDE) then @@ -48,7 +49,10 @@ double precision function comp_source_time_function(t,hdur,it_index) else if (EXTERNAL_SOURCE_TIME_FUNCTION) then ! external stf comp_source_time_function = comp_source_time_function_ext(it_index) - else + elseif(USE_SINSQ_STF)then + ! sin squared stf + comp_source_time_function = comp_source_time_function_ssq(t,hdur) + else ! quasi Heaviside comp_source_time_function = comp_source_time_function_heavi(t,hdur) endif @@ -77,6 +81,39 @@ double precision function comp_source_time_function_heavi(t,hdur) end function comp_source_time_function_heavi +! +!------------------------------------------------------------------------------------------------- +! + + double precision function comp_source_time_function_ssq(t,hdur) + use constants, only: PI + + implicit none + + double precision,intent(in) :: t,hdur + + ! Squared sinusoid STF used in PEGS. Avoids Gaussian in which tiny bits of energy are released + ! from -1.5t0. Energy release only starts at -t0 + ! See, for example, https://doi.org/10.1016/j.epsl.2020.116150 where it is defined as + ! dM/dt = M0/tau sin^2(pi t/(2tau)) where tau is the hdur (not hdur gaussian) + ! Note that this starts from -hdur instead of 0 as defined in attached DOI + + + if(t > -hdur .and. t < hdur)then + comp_source_time_function_ssq = 0.5d0 + ((1/(2.0d0*PI*hdur))*( (hdur * sin(PI*t/hdur)) + PI*t)) + else + ! For regions outside the sinusoid: + if(t.le.-hdur)then + comp_source_time_function_ssq = 0.0d0 + else + comp_source_time_function_ssq = 1.0d0 + endif + endif + + end function comp_source_time_function_ssq + + +! ! !------------------------------------------------------------------------------------------------- diff --git a/src/specfem3D/compute_add_sources.f90 b/src/specfem3D/compute_add_sources.f90 index 9096e3838..4db9c3316 100644 --- a/src/specfem3D/compute_add_sources.f90 +++ b/src/specfem3D/compute_add_sources.f90 @@ -491,7 +491,7 @@ double precision function get_stf_viscoelastic(time_source_dble,isource,it_index ! returns source time function value for specified time - use specfem_par, only: USE_FORCE_POINT_SOURCE,USE_MONOCHROMATIC_CMT_SOURCE,force_stf,hdur,hdur_Gaussian + use specfem_par, only: USE_FORCE_POINT_SOURCE,USE_MONOCHROMATIC_CMT_SOURCE,force_stf,hdur,hdur_Gaussian,USE_SINSQ_STF implicit none @@ -539,7 +539,9 @@ double precision function get_stf_viscoelastic(time_source_dble,isource,it_index if (USE_MONOCHROMATIC_CMT_SOURCE) then f0 = 1.d0 / hdur(isource) ! using half duration as a PERIOD just to avoid changing CMTSOLUTION file format stf = comp_source_time_function_mono(time_source_dble,f0) - else + elseif (USE_SINSQ_STF) then + stf = comp_source_time_function(time_source_dble,hdur(isource),it_index) + else stf = comp_source_time_function(time_source_dble,hdur_Gaussian(isource),it_index) endif endif diff --git a/src/specfem3D/compute_seismograms.F90 b/src/specfem3D/compute_seismograms.F90 index e04a96041..614d29ee9 100644 --- a/src/specfem3D/compute_seismograms.F90 +++ b/src/specfem3D/compute_seismograms.F90 @@ -26,14 +26,14 @@ !===================================================================== - subroutine compute_seismograms(nglob,displ,seismo_current,seismograms) + subroutine compute_seismograms(nglob,displ,seismo_current,seismograms, scale) use constants_solver use specfem_par, only: & nlength_seismogram, & nrec_local,nu_rec,ispec_selected_rec,number_receiver_global, & - scale_displ,hxir_store,hetar_store,hgammar_store + hxir_store,hetar_store,hgammar_store use specfem_par_crustmantle, only: ibool_crust_mantle @@ -41,7 +41,7 @@ subroutine compute_seismograms(nglob,displ,seismo_current,seismograms) integer,intent(in) :: nglob real(kind=CUSTOM_REAL), dimension(NDIM,nglob),intent(in) :: displ - + double precision, intent(in) :: scale integer,intent(in) :: seismo_current real(kind=CUSTOM_REAL), dimension(NDIM,nrec_local,nlength_seismogram),intent(out) :: & @@ -81,9 +81,9 @@ subroutine compute_seismograms(nglob,displ,seismo_current,seismograms) ! store North, East and Vertical components ! distinguish between single and double precision for reals - seismograms(:,irec_local,seismo_current) = real(scale_displ*(nu_rec(:,1,irec)*uxd + & - nu_rec(:,2,irec)*uyd + & - nu_rec(:,3,irec)*uzd), & + seismograms(:,irec_local,seismo_current) = real(scale*(nu_rec(:,1,irec)*uxd + & + nu_rec(:,2,irec)*uyd + & + nu_rec(:,3,irec)*uzd), & kind=CUSTOM_REAL) enddo diff --git a/src/specfem3D/rules.mk b/src/specfem3D/rules.mk index 40461bd97..a6bf77be2 100644 --- a/src/specfem3D/rules.mk +++ b/src/specfem3D/rules.mk @@ -454,6 +454,7 @@ $O/SIEM_solve.solverstatic.o: $O/SIEM_math_library.shared.o $O/SIEM_poisson.solv $O/SIEM_solver_mpi.solverstatic.o: $O/SIEM_math_library.shared.o $O/SIEM_solver_petsc.solverstatic.o: $O/SIEM_math_library.shared.o $O/SIEM_compute_seismograms.solverstatic.o: $O/SIEM_math_library.shared.o +$O/prepare_gravity.solverstatic.o: $O/SIEM_math_library.shared.o ### diff --git a/src/specfem3D/setup_sources_receivers.f90 b/src/specfem3D/setup_sources_receivers.f90 index 307bb3e8d..ecacdc756 100644 --- a/src/specfem3D/setup_sources_receivers.f90 +++ b/src/specfem3D/setup_sources_receivers.f90 @@ -816,6 +816,15 @@ subroutine setup_stf_constants() endif endif + ! Sine-squared STF + if (USE_SINSQ_STF) then + t0 = 0.d0 + if(NSOURCES > 1) write(*,*)'WARNING: USE_SINSQ_STF not tested for NSOURCES > 1' + do isource = 1,NSOURCES + t0 = - min(t0,tshift_src(isource) - hdur(isource)) + enddo + endif + ! checks if user set USER_T0 to fix simulation start time ! note: USER_T0 has to be positive if (USER_T0 > 0.d0) then @@ -1740,6 +1749,12 @@ subroutine setup_receivers_precompute_intp() if (SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3) then allocate(seismograms(NDIM,nrec_local,nlength_seismogram),stat=ier) if (ier /= 0) stop 'Error while allocating seismograms' + + if (FULL_GRAVITY_VAL) then + allocate(seismograms_a(NDIM,nrec_local,nlength_seismogram),stat=ier) + if (ier /= 0) stop 'Error while allocating seismograms_a' + endif + else ! adjoint seismograms allocate(seismograms(NDIM*NDIM,nrec_local,nlength_seismogram),stat=ier) @@ -1760,6 +1775,9 @@ subroutine setup_receivers_precompute_intp() ! initializes seismograms seismograms(:,:,:) = 0._CUSTOM_REAL + if (FULL_GRAVITY_VAL) then + seismograms_a(:,:,:) = 0._CUSTOM_REAL + endif else ! dummy arrays diff --git a/src/specfem3D/specfem3D_par.F90 b/src/specfem3D/specfem3D_par.F90 index c034fe3cc..76c39438d 100644 --- a/src/specfem3D/specfem3D_par.F90 +++ b/src/specfem3D/specfem3D_par.F90 @@ -378,6 +378,8 @@ module specfem_par ! seismograms real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: seismograms + ! acceleration seismograms + real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: seismograms_a integer :: nlength_seismogram integer :: seismo_offset, seismo_current @@ -1136,7 +1138,7 @@ module specfem_par_full_gravity implicit none ! non-dimensionalization - double precision :: scale_accel,scale_pgrav + double precision :: scale_accel, scale_pgrav, scale_phi, scale_hgrav ! seismograms ! perturbed gravitational potential @@ -1385,6 +1387,10 @@ module specfem_par_full_gravity ! for the adjoint Euler scheme for rotation in linear indexing for 3-GLLX points real(kind=CUSTOM_REAL),dimension(:,:), allocatable :: b_A_array_rotationL3, b_B_array_rotationL3 + + + ! testing + logical :: output_background_grav_siem_for_testing = .false. end module specfem_par_full_gravity diff --git a/src/specfem3D/write_seismograms.f90 b/src/specfem3D/write_seismograms.f90 index fbac05cc1..f24bdb553 100644 --- a/src/specfem3D/write_seismograms.f90 +++ b/src/specfem3D/write_seismograms.f90 @@ -35,7 +35,7 @@ subroutine write_seismograms() nrec_local,number_receiver_global,ispec_selected_rec,ispec_selected_source, & it,it_end, & seismo_current,seismo_offset, & - seismograms, & + seismograms, seismograms_a, & nlength_seismogram, & NTSTEP_BETWEEN_OUTPUT_SEISMOS,NTSTEP_BETWEEN_OUTPUT_SAMPLE, & do_save_seismograms, & @@ -44,11 +44,13 @@ subroutine write_seismograms() moment_der,sloc_der,shdur_der,stshift_der, & scale_displ - use specfem_par_crustmantle, only: displ_crust_mantle,b_displ_crust_mantle, & + use specfem_par_crustmantle, only: displ_crust_mantle,accel_crust_mantle,b_displ_crust_mantle, & eps_trace_over_3_crust_mantle,epsilondev_xx_crust_mantle,epsilondev_xy_crust_mantle,epsilondev_xz_crust_mantle, & epsilondev_yy_crust_mantle,epsilondev_yz_crust_mantle, & ibool_crust_mantle + use specfem_par_full_gravity, only: scale_accel + implicit none ! HDF5 is not implemented because ASDF is already available @@ -103,8 +105,12 @@ subroutine write_seismograms() case (1) ! forward run if (.not. GPU_MODE) then - ! on CPU - call compute_seismograms(NGLOB_CRUST_MANTLE,displ_crust_mantle,seismo_current,seismograms) + ! on CPU + call compute_seismograms(NGLOB_CRUST_MANTLE,displ_crust_mantle,seismo_current,seismograms,scale_displ) + ! Full gravity useful to have acceleration outputted + if (FULL_GRAVITY_VAL) then + call compute_seismograms(NGLOB_CRUST_MANTLE,accel_crust_mantle,seismo_current,seismograms_a,scale_accel) + endif else ! on GPU call compute_seismograms_gpu(Mesh_pointer,seismograms,seismo_current,it,it_end,scale_displ,nlength_seismogram) @@ -123,10 +129,10 @@ subroutine write_seismograms() ! on CPU if (OUTPUT_ADJOINT_WAVEFIELD_SEISMOGRAMS) then ! uncomment to output adjoint wavefield instead for seismogram output - call compute_seismograms(NGLOB_CRUST_MANTLE_ADJOINT,displ_crust_mantle,seismo_current,seismograms) + call compute_seismograms(NGLOB_CRUST_MANTLE_ADJOINT,displ_crust_mantle,seismo_current,seismograms, scale_displ) else ! default, backward reconstructed wavefield seismos - call compute_seismograms(NGLOB_CRUST_MANTLE_ADJOINT,b_displ_crust_mantle,seismo_current,seismograms) + call compute_seismograms(NGLOB_CRUST_MANTLE_ADJOINT,b_displ_crust_mantle,seismo_current,seismograms,scale_displ) endif else ! on GPU @@ -189,6 +195,7 @@ subroutine write_seismograms() ! displacement (seismograms) call write_seismograms_to_file(1) if (FULL_GRAVITY_VAL) then + call write_seismograms_to_file(3) ! ground acceleration call write_seismograms_to_file(5) ! seismograms_phi call write_seismograms_to_file(6) ! seismograms_pgrav call write_seismograms_to_file(7) ! seismograms_grav @@ -283,7 +290,7 @@ subroutine write_seismograms_to_file(istore) case (2) component = 'v' ! velocity - not used yet... case (3) - component = 'a' ! acceleration - not used yet... + component = 'a' ! acceleration case (4) component = 'p' ! pressure - not used yet... case (5) @@ -571,7 +578,7 @@ subroutine write_seismograms_to_file(istore) subroutine get_single_trace(istore,irec_local,one_seismogram) use constants, only: NDIM,CUSTOM_REAL - use specfem_par, only: seismograms,nlength_seismogram,seismo_current + use specfem_par, only: seismograms,seismograms_a, nlength_seismogram,seismo_current use specfem_par_full_gravity, only: seismograms_phi,seismograms_pgrav,seismograms_grav,seismograms_corio implicit none @@ -601,11 +608,9 @@ subroutine get_single_trace(istore,irec_local,one_seismogram) continue case (3) ! acceleration - !do i = 1,seismo_current - ! one_seismogram(:,i) = seismograms_a(:,irec_local,i) - !enddo - ! not used yet... - continue + do i = 1,seismo_current + one_seismogram(:,i) = seismograms_a(:,irec_local,i) + enddo case (4) ! pressure !do i = 1,seismo_current @@ -754,8 +759,17 @@ subroutine write_one_seismogram(one_seismogram,irec,irec_local,ftype_flag,compon ! single component only for pressure & gravitational potential if (istore == 4 .or. istore == 5) then chn = bic(1:2)//'P' + elseif (istore == 2) then + ! Velocity + chn = chn(1:3)//'V' + elseif (istore == 3) then + ! Acceleration + chn = chn(1:3)//'A' endif + + + ! backazimuth rotation if (iorientation == 4 .or. iorientation == 5) then ! calculate backazimuth needed to rotate East and North @@ -822,6 +836,17 @@ subroutine write_one_seismogram(one_seismogram,irec,irec_local,ftype_flag,compon write(sisname_big_file,"(a,'.',a,'.',a3,'.sem')") network_name(irec)(1:length_network_name), & station_name(irec)(1:length_station_name),chn + if(istore == 3) then + + write(sisname,"('/',a,'.',a,'.',a4,'.sem')") network_name(irec)(1:length_network_name), & + station_name(irec)(1:length_station_name),chn + + ! create this name also for the text line added to the unique big seismogram file + write(sisname_big_file,"(a,'.',a,'.',a4,'.sem')") network_name(irec)(1:length_network_name), & + station_name(irec)(1:length_station_name),chn + endif + + ! full gravity seismos add an additional component indicator to the name if (istore > 4) then ! using format: **net**.**sta**.channel.C.GRAV.sem.ascii diff --git a/tests/siem/0.configure_siem.default_make.sh b/tests/siem/0.configure_siem.default_make.sh new file mode 100755 index 000000000..c051cca08 --- /dev/null +++ b/tests/siem/0.configure_siem.default_make.sh @@ -0,0 +1,133 @@ +#!/bin/bash +################################################### + +# executable +var=xspecfem3D + +export LD_LIBRARY_PATH=$PETSC_LIB/:$LD_LIBRARY_PATH +OMPI_MCA_btl_openib_allow_ib=1 + +# configuration parameters +CONF_PARAM="--enable-openmp --with-petsc PETSC_INC=$PETSC_INC PETSC_LIB=$PETSC_LIB" + +################################################### +testdir=`pwd` +me=`basename "$0"` + +#checks if ROOT valid +if [ -z "${ROOT}" ]; then export ROOT=../../ ; fi + +# sets source directory +cd $ROOT/ +srcdir=`pwd` + +cd $testdir/ + +# title +echo >> $testdir/results.log +echo "$me in: $testdir" >> $testdir/results.log +echo >> $testdir/results.log + +#cleanup +rm -rf config.log config.status +rm -rf ./bin ./obj ./setup ./OUTPUT_FILES ./DATA/Par_file ./DATABASES_MPI + +mkdir DATABASES_MPI + +# configuration +# (out-of-source compilation) +echo "configuration: $srcdir/configure ${CONF_PARAM}" >> $testdir/results.log +$srcdir/configure ${CONF_PARAM} >> $testdir/results.log 2>&1 + +# checks exit code +if [[ $? -ne 0 ]]; then + echo >> $testdir/results.log + echo "configuration failed, please check..." >> $testdir/results.log + exit 1 +fi + +# we need to output to console output, otherwise tests will fail by timeout in travis +sed -i "s:IMAIN .*:IMAIN = ISTANDARD_OUTPUT:" setup/constants.h >> $testdir/results.log + +# after configuration, we should have a default Par_file in DATA/ +# checks exit code +if [ ! -e DATA/Par_file ]; then + echo "Error: DATA/Par_file not found" >> $testdir/results.log + echo "configuration failed, please check..." >> $testdir/results.log + exit 1 +fi + +# setup model for tests +# overimposes a fixed model setup for testing +# (check also with test_save.f90 in ../meshfem3D test directory) +sed -i "s:^NCHUNKS .*:NCHUNKS = 6:" DATA/Par_file +sed -i "s:^NPROC_XI .*:NPROC_XI = 2:" DATA/Par_file +sed -i "s:^NPROC_ETA .*:NPROC_ETA = 2:" DATA/Par_file + +sed -i "s:^NEX_XI .*:NEX_XI = 48:" DATA/Par_file +sed -i "s:^NEX_ETA .*:NEX_ETA = 48:" DATA/Par_file + +sed -i "s:^MODEL .*:MODEL = transversely_isotropic_prem_plus_3D_crust_1.0:" DATA/Par_file + +sed -i "s:^ANGULAR_WIDTH_XI_IN_DEGREES .*:ANGULAR_WIDTH_XI_IN_DEGREES = 90.d0:" DATA/Par_file +sed -i "s:^ANGULAR_WIDTH_ETA_IN_DEGREES .*:ANGULAR_WIDTH_ETA_IN_DEGREES = 90.d0:" DATA/Par_file +sed -i "s:^CENTER_LATITUDE_IN_DEGREES .*:CENTER_LATITUDE_IN_DEGREES = 90.d0:" DATA/Par_file +sed -i "s:^CENTER_LONGITUDE_IN_DEGREES .*:CENTER_LONGITUDE_IN_DEGREES = 0.d0:" DATA/Par_file +sed -i "s:^GAMMA_ROTATION_AZIMUTH .*:GAMMA_ROTATION_AZIMUTH = 0.d0:" DATA/Par_file + +sed -i "s:^OCEANS .*:OCEANS = .false.:" DATA/Par_file +sed -i "s:^TOPOGRAPHY .*:TOPOGRAPHY = .false.:" DATA/Par_file +sed -i "s:^ELLIPTICITY .*:ELLIPTICITY = .false.:" DATA/Par_file +sed -i "s:^ATTENUATION .*:ATTENUATION = .false.:" DATA/Par_file +sed -i "s:^GRAVITY .*:GRAVITY = .true.:" DATA/Par_file +sed -i "s:^ROTATION .*:ROTATION = .true.:" DATA/Par_file + +sed -i "s:^FULL_GRAVITY .*:FULL_GRAVITY = .true.:" DATA/Par_file + +sed -i "s:^POISSON_SOLVER .*:POISSON_SOLVER = 1:" DATA/Par_file + + +sed -i "s:ABSORBING_CONDITIONS .*:ABSORBING_CONDITIONS = .false.:" DATA/Par_file + +sed -i "s:USE_SINSQ_STF .*:USE_SINSQ_STF = .true.:" DATA/Par_file +sed -i "s:PRINT_SOURCE_TIME_FUNCTION .*:PRINT_SOURCE_TIME_FUNCTION = .true.:" DATA/Par_file + + + +# checks exit code +if [[ $? -ne 0 ]]; then + echo >> $testdir/results.log + echo "compilation failed, please check..." >> $testdir/results.log + exit 1 +fi + +echo "" >> $testdir/results.log +echo "successful configuration" >> $testdir/results.log + +# single compilation +echo "compilation: $var" >> $testdir/results.log +make clean >> $testdir/results.log 2>&1 +make -j 4 $var >> $testdir/results.log 2>&1 + +# checks exit code +if [[ $? -ne 0 ]]; then + echo >> $testdir/results.log + echo "compilation failed, please check..." >> $testdir/results.log + exit 1 +fi + +# checks binary +if [ ! -e bin/$var ]; then + echo "compilation of $var failed, please check..." >> $testdir/results.log + exit 1 +else + echo "binary exists: $var" >> $testdir/results.log +fi +echo "" >> $testdir/results.log + +#cleanup +rm -rf ./bin/* + +echo "successful compilation" >> $testdir/results.log + + diff --git a/tests/siem/1.test_siem_mesh.sh b/tests/siem/1.test_siem_mesh.sh new file mode 100644 index 000000000..30928501e --- /dev/null +++ b/tests/siem/1.test_siem_mesh.sh @@ -0,0 +1,60 @@ +#!/bin/bash +testdir=`pwd` + +# executable +var=test_siem_mesh + +# title +echo >> $testdir/results.log +echo "test: $var" >> $testdir/results.log +echo >> $testdir/results.log + +echo "directory: `pwd`" >> $testdir/results.log + +# clean +mkdir -p bin +rm -f ./bin/$var + +# single compilation +echo "compilation: $var" >> $testdir/results.log + +make -j8 xmeshfem3D >> $testdir/results.log 2>&1 + +echo "" >> $testdir/results.log + +# check +if [ ! -e ./bin/xmeshfem3D ]; then + echo "compilation of $var failed, please check..." >> $testdir/results.log + exit 1 +else + echo "WARNING: TEST IS SIMPLY BUILDING MESHFEM AT THE MOMENT" +fi + +# runs test + +if [ ! -d $testdir/DATABASES_MPI ]; then + mkdir DATABASES_MPI +fi + + +echo "run: `date`" >> $testdir/results.log +mpirun -np 24 ./bin/xmeshfem3D >> $testdir/results.log 2>$testdir/error.log + +# checks exit code +if [[ $? -ne 0 ]]; then + echo "test failed"; echo "error log:"; cat $testdir/error.log; echo "" + exit 1 +fi + +# checks error output (note: fortran stop returns with a zero-exit code) +if [[ -s $testdir/error.log ]]; then + echo "returned ERROR output:" >> $testdir/results.log + cat $testdir/error.log >> $testdir/results.log + exit 1 +fi +rm -f $testdir/error.log + +#cleanup +rm -f bin/$var +# done +echo "successfully tested: `date`" >> $testdir/results.log diff --git a/tests/siem/2.test_gindex.sh b/tests/siem/2.test_gindex.sh new file mode 100644 index 000000000..69d762f0a --- /dev/null +++ b/tests/siem/2.test_gindex.sh @@ -0,0 +1,6 @@ +#!/bin/bash +# Currently not being tested but we do need it to run + +make -j8 gindex3D + +./bin/xgindex3D 24 \ No newline at end of file diff --git a/tests/siem/5.test_sinsq_stf.sh b/tests/siem/5.test_sinsq_stf.sh new file mode 100644 index 000000000..cdd575d3b --- /dev/null +++ b/tests/siem/5.test_sinsq_stf.sh @@ -0,0 +1,66 @@ +#!/bin/bash +testdir=`pwd` + +# executable +var=test_sinsq_stf + +# title +echo >> $testdir/results.log +echo "test: $var" >> $testdir/results.log +echo >> $testdir/results.log + +echo "directory: `pwd`" >> $testdir/results.log + +# clean +mkdir -p bin +rm -f ./bin/$var + +# single compilation +echo "compilation: $var" >> $testdir/results.log + +make -f $var.makefile $var >> $testdir/results.log 2>&1 + +echo "" >> $testdir/results.log + +# check +if [ ! -e ./bin/$var ]; then + echo "compilation of $var failed, please check..." >> $testdir/results.log + exit 1 +fi + +# checks if DATABASES_MPI files were generated from previous 1.test_siem_mesh.sh test +if [ ! -e ./DATABASES_MPI/proc000000_reg3_solver_data.bin ]; then + echo "files in ./DATABASES_MPI/ folder were not generated yet, please check with previous test 1.test_siem_mesh.sh test_save ..." >> $testdir/results.log + exit 1 +fi + +# checks if DATABASES_MPI files were generated from previous 1.test_siem_mesh.sh test +if [ ! -e ./DATABASES_MPI/gdof1_proc23 ]; then + echo "gdof files in ./DATABASES_MPI/ folder were not generated yet, have you run gindex3D? ..." >> $testdir/results.log + exit 1 +fi + + + +# runs test +echo "run: `date`" >> $testdir/results.log +mpirun -np 24 ./bin/$var >> $testdir/results.log 2>$testdir/error.log + +# checks exit code +if [[ $? -ne 0 ]]; then + echo "test failed"; echo "error log:"; cat $testdir/error.log; echo "" + exit 1 +fi + +# checks error output (note: fortran stop returns with a zero-exit code) +if [[ -s $testdir/error.log ]]; then + echo "returned ERROR output:" >> $testdir/results.log + cat $testdir/error.log >> $testdir/results.log + exit 1 +fi +rm -f $testdir/error.log + +#cleanup +rm -f bin/$var +# done +echo "successfully tested: `date`" >> $testdir/results.log diff --git a/tests/siem/test_siem_mesh.makefile b/tests/siem/test_siem_mesh.makefile new file mode 100644 index 000000000..0d8771e34 --- /dev/null +++ b/tests/siem/test_siem_mesh.makefile @@ -0,0 +1,27 @@ +# includes default Makefile from previous configuration +include Makefile + +# test target +default: test_siem_mesh + +## compilation directories +O := ./obj + +OBJECTS = \ + $O/meshfem3D_par.check_module.o \ + $O/auto_ner.shared.o \ + $O/count_elements.shared.o \ + $O/count_points.shared.o \ + $O/define_all_layers.shared.o \ + $O/get_model_parameters.shared.o \ + $O/get_timestep_and_layers.shared.o \ + $O/param_reader.cc.o \ + $O/read_compute_parameters.shared.o \ + $O/read_parameter_file.shared.o \ + $O/read_value_parameters.shared.o \ + $O/shared_par.shared_module.o \ + $(EMPTY_MACRO) + +test_siem_mesh: + ${FCCOMPILE_CHECK} ${FCFLAGS_f90} -o ./bin/test_siem_mesh test_siem_mesh.f90 -I./obj $(OBJECTS) + diff --git a/tests/siem/test_sinsq_stf.f90 b/tests/siem/test_sinsq_stf.f90 new file mode 100644 index 000000000..b93a48f45 --- /dev/null +++ b/tests/siem/test_sinsq_stf.f90 @@ -0,0 +1,71 @@ +program test_sinsq_stf + + use specfem_par + use specfem_par_movie + use manager_adios + use constants, only: PI + + implicit none + integer :: ier,i + character(len=MAX_STRING_LEN) :: junk + double precision :: time, stf, absval, realval, lochdur + + ! initialize the MPI communicator and start the NPROCTOT MPI processes. + call init_mpi() + call world_rank(myrank) + + if (myrank == 0) print *,'program: test_sinsq_stf' + + ! initializes simulation parameters + call initialize_simulation() + + ! sets up reference element GLL points/weights/derivatives + call setup_GLL_points() + + ! starts reading the databases + call read_mesh_databases() + + ! reads topography & bathymetry & ellipticity + call read_topography_bathymetry() + + ! prepares sources and receivers + call setup_sources_receivers() + + ! Add in test here: + + if (myrank == 0)then + + + open(unit = IIN,file = trim(OUTPUT_FILES)//'/plot_source_time_function.txt', status = 'old', iostat=ier, form='formatted') + if (ier /= 0 ) call exit_mpi(0,'Error opening plot_source_time_function file') + + read(IIN, '(A)')junk + read(IIN, '(A)')junk + read(IIN, '(A)')junk + read(IIN, '(A)')junk + + do i = 1, NSTEP + read(IIN,*)time, stf, absval + + ! Assuming HDUR is 20 seconds + lochdur = 20.0d0 + if (time < lochdur) then + realval = 0.5d0 + ((1/(2.0d0*PI*lochdur))*( (lochdur * sin(PI*time/lochdur)) + PI*time)) + else + realval = 1.0d0 + endif + + if (abs(realval - stf) > 1.e-6) then + print *,'ERROR: stf value is incorrect' + stop 1 + endif + + enddo + endif + + ! done + if (myrank == 0) print *,'test_sinsq_stf done successfully' + + ! stop all the MPI processes, and exit + call finalize_mpi() +end program test_sinsq_stf \ No newline at end of file diff --git a/tests/siem/test_sinsq_stf.makefile b/tests/siem/test_sinsq_stf.makefile new file mode 100644 index 000000000..dfa892ec6 --- /dev/null +++ b/tests/siem/test_sinsq_stf.makefile @@ -0,0 +1,17 @@ +# includes default Makefile from previous configuration +include Makefile + +# test target +default: test_sinsq_stf + +## compilation directories +O := ./obj + +OBJECTS = \ + $(specfem3D_SOLVER_OBJECTS) \ + $(specfem3D_SHARED_OBJECTS) \ + $(EMPTY_MACRO) + +test_sinsq_stf: + ${MPIFCCOMPILE_CHECK} ${FCFLAGS_f90} -o ./bin/test_sinsq_stf test_sinsq_stf.f90 -I./obj $(OBJECTS) $(MPILIBS) +