[CIG-SEISMO] [SPECFEM3D_GLOBE] Instability at high(er) resolutions

[Elliott Sales de Andrade]

Hello,

I have been running at a somewhat higher resolution of NEX=1088, as in
the attached Par_file. Unfortunately, the solver blows up in the solid
after a few hundred timesteps. I know people have reported a similar
error in the summer, but that involved the *fluid* and was supposed to
be fixed.

I have tried several older SVN versions (back to about 2008, I think),
but without any luck. I also tried compiling with no optimizations
whatsoever to see if it was a compiler problem (I'm using ifort 12.1.3).

So far, the only working option I found was to reduce the timestep. In
the `rcp_set_timestep_and_layers` subroutine, I have changed the 5%
"safety margin" up to 50%, which seems to work. This reduced the time
step from 0.0475 to 0.025.

I am not certain as to the origin of the DT values in that subroutine,
but I assume they are just pre-computed from some theoretical equation.
While perusing the source, I did notice an additional subroutine in a
different file, `auto_time_stepping`, but I don't believe it is called
for the Par_file I am using.

However, I did notice that the auto_time_stepping subroutine uses the
constant P_VELOCITY_MAX = 11.02827d0. Perhaps this value would be true
for PREM, but it's not true for AK135, and possibly other models. If
this value in auto_time_stepping is incorrect, perhaps the pre-computed
DT values in rcp_set_timestep_and_layers are also incorrect?

-- 
Elliott Sales de Andrade <esalesde at physics.utoronto.ca>
University of Toronto (Physics)

-------------- next part --------------

# forward or adjoint simulation
SIMULATION_TYPE                 = 1
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                         = 2

# angular width of the first chunk (not used if full sphere with six chunks)
ANGULAR_WIDTH_XI_IN_DEGREES   = 90.d0      # angular size of a chunk
ANGULAR_WIDTH_ETA_IN_DEGREES  = 90.d0
CENTER_LATITUDE_IN_DEGREES    = 43.d0
CENTER_LONGITUDE_IN_DEGREES   = -93.d0
GAMMA_ROTATION_AZIMUTH        = 106.1716d0

# 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                          = 1088
NEX_ETA                         = 1088

# number of MPI processors along the two sides of the first chunk
NPROC_XI                        = 34
NPROC_ETA                       = 34

# 1D models with real structure:
# 1D_isotropic_prem, 1D_transversely_isotropic_prem, 1D_iasp91, 1D_1066a, 1D_ak135, 1D_ref, 1D_ref_iso, 1D_jp3d,1D_sea99
#
# 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_ak135_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, s29ea,sea99_jp3d1994,sea99,jp3d1994,heterogen
#
# 3D models with 1D crust: append "_1Dcrust" the 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                           = 3D_anisotropic

# parameters describing the Earth model
OCEANS                          = .false.
ELLIPTICITY                     = .false.
TOPOGRAPHY                      = .false.
GRAVITY                         = .false.
ROTATION                        = .false.
ATTENUATION                     = .true.

# absorbing boundary conditions for a regional simulation
ABSORBING_CONDITIONS            = .true.

# record length in minutes
RECORD_LENGTH_IN_MINUTES        = 30.0d0

# 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                    = .false.
NTSTEP_BETWEEN_FRAMES           = 100
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, 2 or 3, 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

# interval at which we output time step info and max of norm of displacement
NTSTEP_BETWEEN_OUTPUT_INFO      = 1000

# interval in time steps for temporary writing of seismograms
NTSTEP_BETWEEN_OUTPUT_SEISMOS   = 5000
NTSTEP_BETWEEN_READ_ADJSRC      = 1000

# 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.

# rotate seismograms to Radial-Transverse-Z or use default North-East-Z reference frame
ROTATE_SEISMOGRAMS_RT           = .false.

# decide if master process writes all the seismograms or if all processes do it in parallel
WRITE_SEISMOGRAMS_BY_MASTER     = .true.

# save all seismograms in one large combined file instead of one file per seismogram
# to avoid overloading shared non-local file systems such as GPFS for instance
SAVE_ALL_SEISMOS_IN_ONE_FILE    = .true.
USE_BINARY_FOR_LARGE_FILE       = .false.

# flag to impose receivers at the surface or allow them to be buried
RECEIVERS_CAN_BE_BURIED         = .true.

# print source time function
PRINT_SOURCE_TIME_FUNCTION      = .false.