SPectral Elements in Elastodynamics with Discontinuous Galerkin
Tutorials

The folder "TUTORIAL.zip" downloadable from here. contains different examples of elastic wave problems solved with SPEED. The meshes for the following examples have been created using the CUBIT software. See Input files generation for the mesh generation of the following test cases.

Plane wave analysis

Physical problem: plane wave propagation

domain_plane_wave.png
Computational domain considered for the plane wave problem
  • Plane wave with x-polarization and unitary amplitude given on MATE 3
  • Time evolution of the plane wave described by a Ricker wavelet, with peak frequency 1Hz
  • Absorbing boundary conditions imposed on the bottom of the domain
  • Homogeneous Dirichelet boundary conditions (in y and z direction) imposed on the lateral surfaces
  • Homogeneous Neumann boundary condition (free surface) imposed on the top of the domain

Test SEm

The folder PLANE_WAVE/1_TEST_SEM contains:

  • test1_conf_sem.mesh: Mesh constitute by 40 hexahedral elements (12 for MATE 1, 24 for MATE 2 and 4 for MATE 3)
    1  1   hex  1  2  3  4  5  6  7  8
    2  1   hex  2  9  10  3  6  11  12  7
    3  1   hex  4  3  13  14  8  7  15  16
    4  1   hex  3  10  17  13  7  12  18  15
    5  1   hex  5  6  7  8  19  20  21  22
    6  1   hex  6  11  12  7  20  23  24  21
    7  1   hex  8  7  15  16  22  21  25  26
    8  1   hex  7  12  18  15  21  24  27  25
    9  1   hex  19  20  21  22  28  29  30  31
    10  1   hex  20  23  24  21  29  32  33  30
    11  1   hex  22  21  25  26  31  30  34  35
    12  1   hex  21  24  27  25  30  33  36  34
    13  2   hex  35  31  30  34  37  38  39  40
    14  2   hex  31  28  29  30  38  41  42  39
    15  2   hex  34  30  33  36  40  39  43  44
    16  2   hex  30  29  32  33  39  42  45  43
    17  2   hex  37  38  39  40  46  47  48  49
    18  2   hex  38  41  42  39  47  50  51  48
    19  2   hex  40  39  43  44  49  48  52  53
    20  2   hex  39  42  45  43  48  51  54  52
    ...
    
    and 84 boundary faces (4 for absorbing boundary conditions and 80 for Dirichelt boundary conditions)
    1  4  quad  98  94  93  97
    2  4  quad  94  91  92  93
    3  4  quad  97  93  96  99
    4  4  quad  93  92  95  96
    5  5  quad  99  96  83  86
    6  5  quad  96  95  82  83
    7  5  quad  86  83  74  77
    8  5  quad  83  82  73  74
    9  5  quad  77  74  65  68
    10  5  quad  74  73  64  65
    ...
    
  • test1_conf_sem.mate: Description of mechanical properties for MATE 1,2 and 3,
           id  degree     rho        vs     vp       qs  qp
    MATE    1     4       1800      300    600       30  60
    MATE    2     4       2200     2000   4000      200 400
    MATE    3     4       2200     2000   4000      200 400
    
    absorbing and Dirichelet boundary conditions
    ABSO 4 
    DIRY 5 8 0.0 0.0 0.0 0.0
    DIRZ 5 8 0.0 0.0 0.0 0.0
    FUNC 8 0 0.0 
    
    and plane wave forcing term
    PLAX 9 3 1.0
    FUNC 9 4 9.8696 2.0 
    
Note
Polynomial approximation degrees for MATE 1,2 and 3 must be the same if you are using Spectral Element approximation
  • LS.input. File containing physical coordinates of monitored point.
    4
    01  50 50 0 
    02  50 50 -150 
    03  50 50 -250 
    04  50 50 -300
  • SPEED.input: header file
  • FILES_MPI: folder containing *.mpi files, created after partitoning.
  • MONITORS: folder containing output files.
  • run_speed.sh: example of script file for the execution of SPEED on your machine (to be modified according to your installation path). To run SPEED open a terminal, go to "1_TEST_SEM" and digit the command
    ./run_speed.sh 

Test SEm non-linear material

The folder PLANE_WAVE/2_TEST_SEM_NLE contains:

  • test1_conf_sem.mesh: see description in Test SEm
  • test1_conf_sem_nle.mate: Description of mechanical properties for MATE 1,2 and 3,
           id  degree     rho        vs     vp       qs  qp
    MATE    1     4       1800      300    600       30  60
    MATE    2     4       2200     2000   4000      200 400
    MATE    3     4       2200     2000   4000      200 400
    
    absorbing and Dirichelet boundary conditions
    ABSO 4 
    DIRY 5 8 0.0 0.0 0.0 0.0
    DIRZ 5 8 0.0 0.0 0.0 0.0
    FUNC 8 0 0.0 
    
    plane wave forcing term
    PLAX 9 3 1.0
    FUNC 9 4 9.8696 2.0 
    
    and non-linear curve for G/G0 and damping
    MATN   1   4   11  50
    FPEK 0.67 
      Modulus curve
    FUNC 11 60 25   0.000001000   1.000000000  0.000001873   0.999999000  0.000003769   0.998883000  0.000007040   0.998883000  0.000010000   0.985000000  0.000022461   0.955000000  0.000040120   0.928000000  0.000068536   0.860000000  0.000093665   0.825000000  0.000100000   0.800000000  0.000172359   0.700000000  0.000206041   0.650000000  0.000273333   0.600000000  0.000362603   0.520000000  0.000488237   0.450000000  0.000638131   0.370000000  0.000797643   0.320000000  0.001000000   0.280000000  0.001139850   0.260000000  0.001403740   0.230000000  0.002000000   0.190000000  0.003000000   0.150000000  0.007000000   0.105000000  0.007500000   0.104000000  0.010000000   0.100000000
      Damping curve
    FUNC 11 61 25   0.000001000   0.015000000  0.000001873   0.015100000  0.000003769   0.015200000  0.000007040   0.015300000  0.000010000   0.015400000  0.000022461   0.015500000  0.000040120   0.018000000  0.000068536   0.024000000  0.000093665   0.027000000  0.000100000   0.029000000  0.000172359   0.037000000  0.000206041   0.040000000  0.000273333   0.047000000  0.000362603   0.055000000  0.000488237   0.067000000  0.000638131   0.075000000  0.000797643   0.090000000  0.001000000   0.100000000  0.001139850   0.105000000  0.001403740   0.120000000  0.002000000   0.140000000  0.003000000   0.160000000  0.007000000   0.195000000  0.007500000   0.200000000  0.010000000   0.220000000  
    
  • SPEED.input, FILES_MPI, MONITORS: see description in Test SEm
  • run_speed.sh: example of script file for the execution of SPEED on your machine (to be modified according to your installation path). To run SPEED open a terminal, go to "2_TEST_SEM_NLE" and digit the command
    ./run_speed.sh 

Test SEm not-honoring

The folder PLANE_WAVE/INPUTS/3_TEST_SEM_NH contains:

  • test1_conf_sem.mesh: see description in Test SEm
  • test1_conf_sem_nle.mate: Description of mechanical properties for MATE 1,2 and 3,
           id  degree     rho        vs     vp       qs  qp
    MATE    1     4       2200     2000   4000      200 400
    MATE    2     4       2200     2000   4000      200 400
    MATE    3     4       2200     2000   4000      200 400
    
    absorbing and Dirichelet boundary conditions
    ABSO 4 
    DIRY 5 8 0.0 0.0 0.0 0.0
    DIRZ 5 8 0.0 0.0 0.0 0.0
    FUNC 8 0 0.0 
    
    plane wave forcing term
    PLAX 9 3 1.0
    FUNC 9 4 9.8696 2.0 
    
    and case description for not-honoring technique
    CASE 99 1 10
    
    Note
    material properties for MATE 1 are equal to those of MATE 2. The not-honoring method will change runtime the properties of MATE 1 according to the CASE 99 and the definitio of ALL.out and XYZ.out surfaces.
  • SPEED.input, FILES_MPI, MONITORS: see description in Test SEm
  • ALL.out: triangular mesh describing the bottom surface from which the not-honoring method is applied
  • XYZ.out: triangular mesh describing the top surface to which the not-honoring method is applied
  • run_speed.sh: example of script file for the execution of SPEED on your machine (to be modified according to your installation path). To run SPEED open a terminal, go to "3_TEST_SEM_NH" and digit the command
    ./run_speed.sh 

Test DG

The folder PLANE_WAVE/4_TEST_DG contains:

  • test1_nnconf_dg.mesh: Mesh constitute by 124 hexahedra elements (96 for MATE 1, 24 for MATE 2 and 4 for MATE 3)
    1  1   hex  1  2  3  4  5  6  7  8
    2  1   hex  2  9  10  3  6  11  12  7
    3  1   hex  9  13  14  10  11  15  16  12
    4  1   hex  13  17  18  14  15  19  20  16
    5  1   hex  4  3  21  22  8  7  23  24
    6  1   hex  3  10  25  21  7  12  26  23
    7  1   hex  10  14  27  25  12  16  28  26
    8  1   hex  14  18  29  27  16  20  30  28
    9  1   hex  22  21  31  32  24  23  33  34
    10  1   hex  21  25  35  31  23  26  36  33
    ...
    
    and 184 boundary faces (4 for absorbing condition, 28 for DG interface conditions and 152 for Dirichelt boundary conditions )
    1  4  quad  254  246  245  253
    2  4  quad  246  243  244  245
    3  4  quad  253  245  250  256
    4  4  quad  245  244  249  250
    5  5  quad  256  250  248  255
    6  5  quad  250  249  247  248
    7  5  quad  234  231  222  225
    8  5  quad  231  230  221  222
    9  5  quad  225  222  213  216
    10  5  quad  222  221  212  213
    ...
    
    Note
    Surface 6 and 7 define the interface between MATE 3 and MATE 2 respectively. Surface 8 and 9 define the interface between MATE 2 and MATE 1 respectively.
  • test1_nnconf_dg.mate: Description of mechanical properties for MATE 1,2 and 3
           id  degree      rho       vs     vp       qs  qp
    MATE    1     4       1800      300    600       30  60
    MATE    2     4       2200     2000   4000      200 400
    MATE    3     4       2200     2000   4000      200 400
    
    absorbing and Dirichelet boundary conditions
    ABSO 4 
    DIRY 5 8 0.0 0.0 0.0 0.0
    DIRZ 5 8 0.0 0.0 0.0 0.0
    FUNC 8 0 0.0 
    
    DG interface conditions
    DGIC 6  1
    DGIC 7  0
    DGIC 8  1
    DGIC 9  0
    
    and plane wave forcing term
    PLAX 9 3 1.0
    FUNC 9 4 9.8696 2.0 
    
    Note
    Polynomial degrees for MATE 1,2 and 3 can be set independently one from the other for Discontinuous Galerkin Spectral Element approxiamtion.
  • LS.input. File containing physical coordinates of monitored point.
    4
    01  50 50 0 
    02  50 50 -150 
    03  50 50 -250 
    04  50 50 -300
    
  • SPEED.input: header file
  • FILES_MPI: folder containing *.mpi files, created after partitoning.
  • MONITORS: folder containing output files.
  • run_speed.sh: example of script file for the execution of SPEED on your machine (to be modified according to your installation path). To run SPEED open a terminal, go inside the "4_TEST_DG" folder and digit the command
    ./run_speed.sh 

Post-processing and visualization

The post-processing and visualization of the results is obtained through MATLAB scripts, provided in the folder POST-PROC inside each test case folder.

  • Use the MATLAB script REWRITE_MONITOR_FORMAT.m to rewrite the output files (from MONITORXXXXXX.D to monitorxxxxxx.d) stored in the directory MONITORS.
  • Run the MATLAB script PLOT_MONITOR.m to plot the output results.
Attention
Modify, if necessary, the part TO BE DEFINED BY USER before running the scripts.
plane_wave.png
Comparison between the 4 different models and methods presented in the planve wave tutorial. Time-hhistory of the x-displacement recorded at point (50,50,0)

.

LOH1 test case

Physical problem: Point source wave propagation

Benchmark tectcase Layer over a half space (LOH1) described in:
"Day S.M., Bradley C.R. 2001. Memory-efficient simulation of anelastic wave propagation. Bulletin of the Seismological Society of America, 91(3): 520–531."

  • The folder TEST_CASE_LOH1/INPUTS contains two subdirectories DGSEM and SEM in which input file are prepared to tackle the problem with Discontinuous or Continuous Spectral elements, respectively.
  • Input files are prepared according to the paper "I.Mazzieri, M.Stupazzini, R.Guidotti, C.Smerzini: SPEED: SPectral Elements in Elastodynamics with Discontinuous Galerkin: a non-conforming approach for 3D multi-scale problems, Int. J. Numer. Meth. Eng., (2013); 95(12): 991-1010."
test_case_loh1.png
LOH1 conforming model, having size of 30 × 30 × 17 km: (Top) conforming mesh with 814.833 hexahedral elements, varying from size of 100 m, in the first quadrant, to 300 m in the remaining part of the domain; (Bottom) non-conforming model with 70.228 hexahedral elements, having size of around 400 m in the upper layer (1 km thickness) and size of around 650 m in the lower layer (16 km thickness).

Execution

  • The directories TEST_CASE_LOH1/INPUTS/DGSEM and TEST_CASE_LOH1/INPUTS/SEM contain an example of input files

Post-processing and visualization

  • Copy all files contained in TEST_CASE_LOH1/INPUTS/DGSEM/MONITORS or TEST_CASE_LOH1/INPUTS/SEM/MONITORS into the folder POST-PROC/MONITORS
  • Run the MATLAB script MAIN.m

Croissant Valley test case

Physical problem: Plane wave propagation

Benchmark tectcase proposed by F.J. Sanchez-Sesma, F. Luzon 1995. Seismic Response of a Three-Dimensional Alluvial Valleys for Incident P, S and Rayleigh Waves. Bulletin of the Seismological Society of America, 85, 890-899. See also:

  • I. Mazzieri 2012. "Non-Conforming High Order Methods for the Elastodynamics Equation", Ph.D. Thesis, Politecnico di Milano.
  • M. Stupazzini 2004. "A Spectral Element Approach for 3D Dynamic Soil Structure Interaction Problems", Ph.D. Thesis, Politecnico di Milano.

The folder CROISSANT_VALLEY contains the following subdirectories:

  • Inputs: input files for different simulation with SEM, DGSEM and SEM with Not-Honoring.
  • Mesh: mesh files for the different discretization approaches
  • SesmaLuzon: semi-analitycal solutions obtained by IBEM
  • SrcMatlab: source file for the post-processing
  • SesmaLuzonVsSpeed.m: main Matlab program for the post-porcessing
croissantvalley.png
The numerical mesh adopted for the test of plane wave load with angle of incidence orthogonal to the free surface. The numerical mesh has 13038 hexahedral elements and 864 957 nodes with spectral degree 4

.

Execution

  • The directories Inputs/1ConformingSEM, Inputs/2ConformingDG and Inputs/3NotHonoringSEM contain examples of input files

Post-processing and visualization

  • Run the MATLAB script SesmaLuzonVsSPeed.m