simulation issue [solved]

[geostar]

I recently tried running a forward simulation using SPECFEM3D_GLOBE, but encountered the following error:

forward simulation became unstable in solid and blew up
Error detected, aborting MPI... proc 0

I used the s362ani model with CRUST1.0 and set up the following chunk parameters:

ANGULAR_WIDTH_XI_IN_DEGREES = 90.d0
ANGULAR_WIDTH_ETA_IN_DEGREES = 90.d0
CENTER_LATITUDE_IN_DEGREES = 20.d0
CENTER_LONGITUDE_IN_DEGREES = 110.d0
GAMMA_ROTATION_AZIMUTH = 60.d0

NEX_XI = 288
NEX_ETA = 288

I was able to fix the issue by reducing the time step (DT). In the file /shared/get_timestep_and_layers.f90, between lines 479-499, I found this section:

if (REGIONAL_MOHO_MESH) then
   ! hard coded number of crustal element layers and time step

   ! checks
   if (NCHUNKS > 1 ) stop 'regional moho mesh: NCHUNKS error in rcp_set_timestep_and_layers'
   if (HONOR_1D_SPHERICAL_MOHO ) return

   ! original values
   !print *,'NER:',NER_CRUST
   !print *,'DT:',DT

   ! enforce 3 element layers
   NER_CRUST = 3

   ! increased stability, empirical
   DT = DT*(1.d0 + 0.5d0)

   if (REGIONAL_MOHO_MESH_EUROPE ) DT = 0.17 ! Europe
   if (REGIONAL_MOHO_MESH_ASIA ) DT = 0.15 ! Asia & Middle East

 endif

When I removed line 494, which modifies the time step for a three-layer regional Moho mesh

DT = DT*(1.d0 + 0.5d0)

, the simulation ran smoothly. I think this issue might need further attention.

[danielpeter]

thanks for the note!

could you send the mesher output file to see how the mesh evaluation looked liked?

in your setup, there are likely some distorted crustal elements affecting the stability. the DT have all been chosen empirically and as big as possible to avoid unnecessary computational costs. unfortunately, that's not working in all possible setups.

did you know that you could add a line in the Par_file to override the emprical DT value? by setting

DT   = 0.015

with a value you think should be stable, you can re-run the simulation without the need of modifying the source code file.

[geostar]

Thanks for letting me know I can overwrite DT by adding it as a parameter in the Par_file. I thought the code automatically set DT based on velocities and the element size to meet the CFL condition, and it worked well before. This time, I tried a larger chunk that covers different part of the Earth, and that's when the issue came up.

For distorted crustal elements, I assume the code will raise a different error, like the Jacobian being negative?

Please see the mesher output below:

"output_mesher.txt" 406L, 12206B 406,1 Bot

---

*** Specfem3D MPI Mesher ***

---

There are 144 MPI processes
Processes are numbered from 0 to 143

There are 288 elements along xi in each chunk
There are 288 elements along eta in each chunk

There are 12 slices along xi in each chunk
There are 12 slices along eta in each chunk
There is a total of 144 slices in each chunk
There are 1 chunks in the global mesh
There is a total of 144 slices in the global mesh

NGLLX = 5
NGLLY = 5
NGLLZ = 5

Shape functions defined by NGNOD = 27 control nodes
Surface shape functions defined by NGNOD2D = 9 control nodes

model: s362ani
incorporating the oceans using equivalent load
incorporating ellipticity
incorporating surface topography
incorporating self-gravitation (Cowling approximation)
incorporating rotation
incorporating attenuation using 3 standard linear solids

incorporating 3-D lateral variations
no heterogeneities in the mantle
incorporating crustal variations
using one layer only in PREM crust
incorporating anisotropy
no inner-core anisotropy
no general mantle anisotropy

Reference radius of the Earth used is 6371.00000000000 km

Central cube is at a radius of 940.000000000000 km
creating global slice addressing

incorporating topography
topography/bathymetry: min/max = -7747 5507

incorporating crustal model: CRUST1.0

additional mesh optimizations

moho:
enforcing regional 3-layer crust
incorporating element stretching for 3-D moho surface with deep moho

internal topography 410/660:
incorporating element stretching for 3-D internal surfaces

---

creating mesh in region 1
this region is the crust and mantle

---

first pass

...allocating arrays

...setting up layers

...creating mesh elements
creating layer 1 out of 10
10.0% current clock (NOT elapsed) time is: 17h 59min 05sec
creating layer 2 out of 10
20.0% current clock (NOT elapsed) time is: 17h 59min 05sec
creating layer 3 out of 10
30.0% current clock (NOT elapsed) time is: 17h 59min 05sec
creating layer 4 out of 10
40.0% current clock (NOT elapsed) time is: 17h 59min 05sec
creating layer 5 out of 10
50.0% current clock (NOT elapsed) time is: 17h 59min 05sec
creating layer 6 out of 10
60.0% current clock (NOT elapsed) time is: 17h 59min 05sec
creating layer 7 out of 10
70.0% current clock (NOT elapsed) time is: 17h 59min 05sec
creating layer 8 out of 10
80.0% current clock (NOT elapsed) time is: 17h 59min 05sec
creating layer 9 out of 10
90.0% current clock (NOT elapsed) time is: 17h 59min 05sec
creating layer 10 out of 10
100.0% current clock (NOT elapsed) time is: 17h 59min 06sec

...creating global addressing

...creating MPI buffers

second pass

...allocating arrays

...setting up layers

...creating mesh elements
creating layer 1 out of 10
10.0% current clock (NOT elapsed) time is: 17h 59min 10sec
creating layer 2 out of 10
20.0% current clock (NOT elapsed) time is: 17h 59min 11sec
creating layer 3 out of 10
30.0% current clock (NOT elapsed) time is: 17h 59min 12sec
creating layer 4 out of 10
40.0% current clock (NOT elapsed) time is: 17h 59min 14sec
creating layer 5 out of 10
50.0% current clock (NOT elapsed) time is: 17h 59min 19sec
creating layer 6 out of 10
60.0% current clock (NOT elapsed) time is: 17h 59min 21sec
creating layer 7 out of 10
70.0% current clock (NOT elapsed) time is: 17h 59min 21sec
creating layer 8 out of 10
80.0% current clock (NOT elapsed) time is: 17h 59min 23sec
creating layer 9 out of 10
90.0% current clock (NOT elapsed) time is: 17h 59min 25sec
creating layer 10 out of 10
100.0% current clock (NOT elapsed) time is: 17h 59min 32sec

...precomputing Jacobian

...creating chunk buffers

----- creating chunk buffers -----

There are 12 slices along xi in each chunk
There are 12 slices along eta in each chunk
There is a total of 144 slices in each chunk
There are 1 chunks
There is a total of 144 slices in all the chunks

There is a total of 12 messages to assemble faces between chunks

only one chunk, no need to create chunk buffers

...preparing MPI interfaces

crust/mantle region:
#max of points in MPI buffers along xi npoin2D_xi = 9285
#max of array elements transferred npoin2D_xi*NDIM = 27855

#max of points in MPI buffers along eta npoin2D_eta = 9285
#max of array elements transferred npoin2D_eta*NDIM = 27855

crust mantle MPI:
maximum interfaces: 8
MPI addressing maximum interfaces: 8
MPI addressing : all interfaces okay

total MPI interface points : 4999764
unique MPI interface points: 4805196
maximum valence : 3
total unique MPI interface points: 4805196

...element inner/outer separation

for overlapping of communications with calculations:

percentage of edge elements in crust/mantle 14.77987 %
percentage of volume elements in crust/mantle 85.22012 %

...element mesh coloring
mesh coloring: F

...creating mass matrix
creates exact mass matrix for rotation
updates mass matrix with Stacey boundary corrections
updates mass matrix with ocean load

...saving binary files

calculated top area:    2.09512003294649

calculated bottom area: 0.625251241238778

...saving AVS or DX mesh files

---

creating mesh in region 2
this region is the outer core

---

first pass

...allocating arrays

...setting up layers

...creating mesh elements
creating layer 1 out of 2
50.0% current clock (NOT elapsed) time is: 17h 59min 40sec
creating layer 2 out of 2
100.0% current clock (NOT elapsed) time is: 17h 59min 40sec

...creating global addressing

...creating MPI buffers

second pass

...allocating arrays

...setting up layers

...creating mesh elements
creating layer 1 out of 2
50.0% current clock (NOT elapsed) time is: 17h 59min 40sec
creating layer 2 out of 2
100.0% current clock (NOT elapsed) time is: 17h 59min 41sec

...precomputing Jacobian

...creating chunk buffers

----- creating chunk buffers -----

There are 12 slices along xi in each chunk
There are 12 slices along eta in each chunk
There is a total of 144 slices in each chunk
There are 1 chunks
There is a total of 144 slices in all the chunks

There is a total of 12 messages to assemble faces between chunks

only one chunk, no need to create chunk buffers

...preparing MPI interfaces

outer core region:
#max of points in MPI buffers along xi npoin2D_xi = 2505
#max of array elements transferred npoin2D_xi*NDIM = 7515

#max of points in MPI buffers along eta npoin2D_eta = 2553
#max of array elements transferred npoin2D_eta*NDIM = 7659

outer core MPI:
maximum interfaces: 8
MPI addressing maximum interfaces: 8
MPI addressing : all interfaces okay

total MPI interface points : 1381620
unique MPI interface points: 1252524
maximum valence : 3
total assembled MPI interface points: 1252524

...element inner/outer separation
percentage of edge elements in outer core 36.52778 %
percentage of volume elements in outer core 63.47222 %

...element mesh coloring
mesh coloring: F

...creating mass matrix
creates exact mass matrix for rotation
updates mass matrix with Stacey boundary corrections

...saving binary files

calculated top area:   0.625251241238778

calculated bottom area: 7.703250977623870E-002

...saving AVS or DX mesh files

---

creating mesh in region 3
this region is the inner core

---

first pass

...allocating arrays

...setting up layers

...creating mesh elements
creating layer 1 out of 1
100.0% current clock (NOT elapsed) time is: 17h 59min 44sec

...creating global addressing

...creating MPI buffers

second pass

...allocating arrays

...setting up layers

...creating mesh elements
creating layer 1 out of 1
100.0% current clock (NOT elapsed) time is: 17h 59min 44sec

...precomputing Jacobian

...creating chunk buffers

----- creating chunk buffers -----

There are 12 slices along xi in each chunk
There are 12 slices along eta in each chunk
There is a total of 144 slices in each chunk
There are 1 chunks
There is a total of 144 slices in all the chunks

There is a total of 12 messages to assemble faces between chunks

only one chunk, no need to create chunk buffers

...preparing MPI interfaces

inner core region:
#max of points in MPI buffers along xi npoin2D_xi = 221
#max of array elements transferred npoin2D_xi*NDIM = 663

#max of points in MPI buffers along eta npoin2D_eta = 221
#max of array elements transferred npoin2D_eta*NDIM = 663

inner core MPI:
maximum interfaces: 8
MPI addressing maximum interfaces: 8
MPI addressing : all interfaces okay

total MPI interface points : 124916
unique MPI interface points: 108460
maximum valence : 3
total assembled MPI interface points: 108460

...element inner/outer separation
percentage of edge elements in inner core 55.55556 %
percentage of volume elements in inner core 44.44444 %
...element mesh coloring
mesh coloring: F

...creating mass matrix
creates exact mass matrix for rotation
updates mass matrix with Stacey boundary corrections

...saving binary files

calculated top area:   7.703250977623870E-002

calculated bottom area: 4.561931621663209E-002

...saving AVS or DX mesh files

calculated volume: 0.696184682988442

Repartition of elements in regions:

total number of elements in each slice: 8388

• crust and mantle: 90.98712 %
• outer core: 8.583691 %
• inner core: 0.4291846 %

for some mesh statistics, see comments in file OUTPUT_FILES/values_from_mesher.
h

Load balancing = 100 % by definition

the time step of the solver will be DT = 0.1325250

using single precision for the calculations

smallest and largest possible floating-point numbers are: 1.1754944E-38
3.4028235E+38

Elapsed time for mesh generation and buffer creation in seconds =
50.0976700587198
Elapsed time for mesh generation and buffer creation in hh:mm:ss = 0 h 00 m 50 s

End of mesh generation

                                                            406,1         Bot

[danielpeter]

thanks for the output!

I notice that you use an older code version. could you also use the most recent devel branch version to test if this stability issue still persists? and if so, please send again the output_mesher.txt file :)