[CIG-SEISMO] Floating-point exception - erroneous arithmetic operation
Moritz Wehler
moritzwehler at muenchen-ist-toll.de
Wed May 27 06:36:37 PDT 2015
Hey there,
I have some issues with increasing the number of elements along the
z-axis with specfem2D.
As soon as i double the number of elements from 111 (in total along the
z-axis) to 222 (in total), I get the following error message:
"Program received signal SIGFPE: Floating-point exception - erroneous
arithmetic operation.
Backtrace for this error:
#0 0x7F62A8D0D7D7
#1 0x7F62A8D0DDDE
#2 0x7F62A8449D3F
#3 0x433DBE in compute_forces_viscoelastic_
#4 0x4DF021 in MAIN__ at specfem2D.F90:?
Gleitkomma-Ausnahme (Speicherabzug geschrieben)"
I should add that everything worked perfectly fine BEFORE I increased
the element number. The corresponding par_file and interface_file are
attached.
Hope I didn't forget anything important and that someone can help me.
Cheers, Mo
-------------- next part --------------
A non-text attachment was scrubbed...
Name: interfaces_simple_topo_curved.dat
Type: application/x-ns-proxy-autoconfig
Size: 965 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/cig-seismo/attachments/20150527/528546f3/attachment-0001.dat>
-------------- next part --------------
# title of job
title = Test for simple topography
# axisymmetric (2.5D) or Cartesian planar (2D) simulation
AXISYM = .false.
# forward or adjoint simulation
SIMULATION_TYPE = 1 # 1 = forward, 3 = adjoint + kernels; 2 is purposely UNUSED (for compatibility with the numbering of our 3D codes)
NOISE_TOMOGRAPHY = 0 # 0 = earthquake simulation, 1/2/3 = noise simulation
SAVE_FORWARD = .false. # save the last frame, needed for adjoint simulation
# parameters concerning partitioning
nproc = 1 # number of processes
partitioning_method = 3 # SCOTCH = 3, ascending order (very bad idea) = 1
ngnod = 9 # number of control nodes per element (4 or 9)
initialfield = .false. # use a plane wave as source or not
add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
assign_external_model = .false. # define external earth model or not
READ_EXTERNAL_SEP_FILE = .false. # Read external model from DATA/model_velocity.dat_input, or use routine
ATTENUATION_VISCOELASTIC_SOLID = .true. # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
ATTENUATION_PORO_FLUID_PART = .true. # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
Q0 = 100 # quality factor for viscous attenuation
freq0 = 10 # frequency for viscous attenuation
p_sv = .true. # set the type of calculation (P-SV or SH/membrane ewaves)
# time step parameters
nt = 2000 # total number of time steps
deltat = 1.d-3 # duration of a time step (see section "How to choose the time step" of the manual for how to do this)
USER_T0 = 0.0d0 # use this t0 as earliest starting time rather than the automatically calculated one
time_stepping_scheme = 1 # 1 = Newmark (2nd order), 2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta), 3 = classical RK4 4th-order 4-stage Runge-Kutta
# acoustic forcing
ACOUSTIC_FORCING = .true. # acoustic forcing of an acoustic medium with a rigid interface
# source parameters
NSOURCES = 1 # number of sources (source info read from DATA/CMTSOLUTION file)
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
# for viscoelastic attenuation
N_SLS = 2 # number of standard linear solids for attenuation
f0_attenuation = 10. # (Hz) relevant only if source is a Dirac or a Heaviside, otherwise it is f0 the dominant frequency of the source in the CMTSOLUTION file
# shift (i.e. change) velocities read from the input file to take average physical dispersion into account,
# i.e. if needed change the reference frequency at which these velocities are defined internally in the code:
# by default, the seismic velocity values that are read at the end of this Par_file of the code are supposed to be the unrelaxed values,
# i.e. the velocities at infinite frequency. We may want to change this and impose that the values read are those for a given frequency (here f0_attenuation).
# (when we do this, velocities will then slightly decrease and waves will thus slightly slow down)
READ_VELOCITIES_AT_f0 = .true.
# receiver set parameters for seismograms
seismotype = 2 # record 1=displ 2=veloc 3=accel 4=pressure 5=curl of displ 6=the fluid potential
NSTEP_BETWEEN_OUTPUT_SEISMOS = 5000000 # every how many time steps we save the seismograms (costly, do not use a very small value; if you use a very large value that is larger than the total number of time steps of the run, the seismograms will automatically be saved once at the end of the run anyway)
save_ASCII_seismograms = .true. # save seismograms in ASCII format or not
save_binary_seismograms_single = .true. # save seismograms in single precision binary format or not (can be used jointly with ASCII above to save both)
save_binary_seismograms_double = .false. # save seismograms in double precision binary format or not (can be used jointly with both flags above to save all)
SU_FORMAT = .false. # output single precision binary seismograms in Seismic Unix format (adjoint traces will be read in the same format)
subsamp_seismos = 1 # subsampling of the seismograms to create smaller files (but less accurately sampled in time)
use_existing_STATIONS = .false. # use an existing STATION file found in ./DATA or create a new one from the receiver positions below in this Par_file
nreceiversets = 4 # number of receiver sets
anglerec = 0.d0 # angle to rotate components at receivers
rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
save_ASCII_kernels = .true. # save sensitivity kernels in ASCII format (much bigger files, but compatible with current GMT scripts) or in binary format
# first receiver set (repeat these 6 lines and adjust nreceiversets accordingly)
nrec = 1 # number of receivers
xdeb = 10289. # first receiver x in meters
zdeb = 9847. # first receiver z in meters
xfin = 3700. # last receiver x in meters (ignored if only one receiver)
zfin = 2200. # last receiver z in meters (ignored if only one receiver)
enreg_surf_same_vertical = .true. # receivers inside the medium or at the surface
# first receiver set (repeat these 6 lines and adjust nreceiversets accordingly)
nrec = 1 # number of receivers
xdeb = 9409. # first receiver x in meters
zdeb = 9863. # first receiver z in meters
xfin = 3700. # last receiver x in meters (ignored if only one receiver)
zfin = 2200. # last receiver z in meters (ignored if only one receiver)
enreg_surf_same_vertical = .true. # receivers inside the medium or at the surface
# first receiver set (repeat these 6 lines and adjust nreceiversets accordingly)
nrec = 1 # number of receivers
xdeb = 8515. # first receiver x in meters
zdeb = 9873. # first receiver z in meters
xfin = 3700. # last receiver x in meters (ignored if only one receiver)
zfin = 2200. # last receiver z in meters (ignored if only one receiver)
enreg_surf_same_vertical = .true. # receivers inside the medium or at the surface
# first receiver set (repeat these 6 lines and adjust nreceiversets accordingly)
nrec = 1 # number of receivers
xdeb = 14879. # first receiver x in meters
zdeb = 9900. # first receiver z in meters
xfin = 3700. # last receiver x in meters (ignored if only one receiver)
zfin = 2200. # last receiver z in meters (ignored if only one receiver)
enreg_surf_same_vertical = .true. # receivers inside the medium or at the surface
# display parameters
NSTEP_BETWEEN_OUTPUT_INFO = 400 # every how many time steps we display information about the simulation (costly, do not use a very small value)
####
NSTEP_BETWEEN_OUTPUT_IMAGES = 400 # every how many time steps we draw JPEG or PostScript pictures of the simulation (costly, do not use a very small value)
cutsnaps = 1. # minimum amplitude kept in % for the JPEG and PostScript snapshots; amplitudes below that are muted
#### for JPEG color images ####
output_color_image = .true. # output JPEG color image of the results every NSTEP_BETWEEN_OUTPUT_IMAGES time steps or not
imagetype_JPEG = 6 # display 1=displ_Ux 2=displ_Uz 3=displ_norm 4=veloc_Vx 5=veloc_Vz 6=veloc_norm 7=accel_Ax 8=accel_Az 9=accel_norm 10=pressure
factor_subsample_image = 1.0d0 # (double precision) factor to subsample color images output by the code (useful for very large models)
USE_CONSTANT_MAX_AMPLITUDE = .true. # by default the code normalizes each image independently to its maximum; use this option to use the global maximum below instead
CONSTANT_MAX_AMPLITUDE_TO_USE = 1.17d-2 # constant maximum amplitude to use for all color images if the above USE_CONSTANT_MAX_AMPLITUDE option is true
POWER_DISPLAY_COLOR = 0.30d0 # non linear display to enhance small amplitudes in JPEG color images
DRAW_SOURCES_AND_RECEIVERS = .true. # display sources as orange crosses and receivers as green squares in JPEG images or not
DRAW_WATER_IN_BLUE = .true. # display acoustic layers as constant blue in JPEG images, because they likely correspond to water in the case of ocean acoustics or in the case of offshore oil industry experiments (if off, display them as greyscale, as for elastic or poroelastic elements, for instance for acoustic-only oil industry models of solid media)
USE_SNAPSHOT_NUMBER_IN_FILENAME = .false. # use snapshot number in the file name of JPEG color snapshots instead of the time step (for instance to create movies in an easier way later)
#### for PostScript snapshots ####
output_postscript_snapshot = .false. # output Postscript snapshot of the results every NSTEP_BETWEEN_OUTPUT_IMAGES time steps or not
imagetype_postscript = 2 # display 1=displ vector 2=veloc vector 3=accel vector; small arrows are displayed for the vectors
meshvect = .false. # display mesh on PostScript plots or not
modelvect = .false. # display velocity model on PostScript plots or not
boundvect = .false. # display boundary conditions on PostScript plots or not
interpol = .false. # interpolation of the PostScript display on a regular grid inside each spectral element, or use the non-evenly spaced GLL points
pointsdisp = 6 # number of points in each direction for interpolation of PostScript snapshots (set to 1 for lower-left corner only)
subsamp_postscript = 1 # subsampling of background velocity model in PostScript snapshots
sizemax_arrows = 1.d0 # maximum size of arrows on PostScript plots in centimeters
US_LETTER = .false. # use US letter or European A4 paper for PostScript plots
#### for wavefield dumps ####
NSTEP_BETWEEN_OUTPUT_WAVE_DUMPS = 400 # every how many time steps we dump results of the simulation as ASCII or binary files (costly, do not use a very small value)
output_wavefield_dumps = .false. # output wave field to a text file every NSTEP_BETWEEN_OUTPUT_TEXT_DUMPS time steps (creates very big files)
imagetype_wavefield_dumps = 2 # display 1=displ vector 2=veloc vector 3=accel vector 4=pressure
use_binary_for_wavefield_dumps = .false. # use ASCII or single-precision binary format for the wave field dumps
####
output_grid_Gnuplot = .false. # generate a GNUPLOT file containing the grid, and a script to plot it
output_grid_ASCII = .false. # dump the grid in an ASCII text file consisting of a set of X,Y,Z points or not
output_energy = .false. # compute and output total acoustic and elastic energy curves (slows down the code significantly)
# velocity and density models
nbmodels = 10 # nb of different models
# define models as
# I: (model_number 1 rho Vp Vs 0 0 QKappa Qmu 0 0 0 0 0 0) or
# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 c12 c23 c25 0 0 0) or
# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
#
# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
#
# For anisotropic media in the most general 2D case we have the following convention for the stress-strain relationship:
# ! implement anisotropy in 2D
# sigma_xx = c11*dux_dx + c13*duz_dz + c15*(duz_dx + dux_dz)
# sigma_zz = c13*dux_dx + c33*duz_dz + c35*(duz_dx + dux_dz)
# sigma_xz = c15*dux_dx + c35*duz_dz + c55*(duz_dx + dux_dz)
# sigma_yy is not equal to zero in the plane strain formulation
# but is used only in post-processing if needed,
# to compute pressure for display or seismogram recording purposes
# sigma_yy = c12*dux_dx + c23*duz_dz + c25*(duz_dx + dux_dz)
# where the notations are for instance duz_dx = d(Uz) / dx.
#
# For anisotropic elastic media the last three parameters, c12 c23 c25, are used only when the user asks the code to compute pressure for display
# or seismogram recording purposes. Thus, if you do not know these parameters for your anisotropic material and/or if you do not plan to display or record pressure you
# can ignore them and set them to zero. When pressure is used these three parameters are needed because the code needs to compute sigma_yy,
# which is not equal to zero in the plane strain formulation.
#(1 rho Vp Vs 0 0 QKappaQmu 0 0 0 0 0 0)
1 1 2500.d0 6000.d0 3640.d0 0 0 9999 9999 0 0 0 0 0 0
2 1 2500.d0 5600.d0 3390.d0 0 0 9999 9999 0 0 0 0 0 0
3 1 2500.d0 4760.d0 2880.d0 0 0 9999 9999 0 0 0 0 0 0
4 1 2500.d0 4340.d0 2630.d0 0 0 9999 9999 0 0 0 0 0 0
5 1 2500.d0 4180.d0 2530.d0 0 0 9999 9999 0 0 0 0 0 0
6 1 2500.d0 3960.d0 2400.d0 0 0 9999 9999 0 0 0 0 0 0
7 1 2500.d0 3710.d0 2250.d0 0 0 9999 9999 0 0 0 0 0 0
8 1 2500.d0 3380.d0 2050.d0 0 0 9999 9999 0 0 0 0 0 0
9 1 2500.d0 2870.d0 1740.d0 0 0 9999 9999 0 0 0 0 0 0
10 1 2500.d0 2180.d0 1320.d0 0 0 9999 9999 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .false.
# absorbing boundary active or not
PML_BOUNDARY_CONDITIONS = .true.
NELEM_PML_THICKNESS = 3
ROTATE_PML_ACTIVATE = .false.
ROTATE_PML_ANGLE = 30.
STACEY_ABSORBING_CONDITIONS = .false.
ADD_SPRING_TO_STACEY = .true.
# for horizontal periodic conditions: detect common points between left and right edges
ADD_PERIODIC_CONDITIONS = .false.
# horizontal periodicity distance for periodic conditions
PERIODIC_HORIZ_DIST = 0.3597d0
#-----------------------------------------------------------------------------
# PARAMETERS FOR EXTERNAL MESHING
# data concerning mesh, when generated using third-party app (more info in README)
# (see also absorbing_conditions above)
mesh_file = ./DATA/Mesh_canyon/canyon_mesh_file # file containing the mesh
nodes_coords_file = ./DATA/Mesh_canyon/canyon_nodes_coords_file # file containing the nodes coordinates
materials_file = ./DATA/Mesh_canyon/canyon_materials_file # file containing the material number for each element
free_surface_file = ./DATA/Mesh_canyon/canyon_free_surface_file # file containing the free surface
axial_elements_file = ./DATA/axial_elements_file # file containing the axial elements if AXISYM is true
absorbing_surface_file = ./DATA/Mesh_canyon/canyon_absorbing_surface_file # file containing the absorbing surface
acoustic_forcing_surface_file = ./DATA/MSH/Surf_acforcing_Bottom_enforcing_mesh # file containing the acoustic forcing surface
CPML_element_file = Elements_CPML_list # file containing the CPML element numbers
tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
#-----------------------------------------------------------------------------
# PARAMETERS FOR INTERNAL MESHING
# file containing interfaces for internal mesh
interfacesfile = ../interfaces_simple_topo_curved.dat
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
xmax = 15000.d0 # abscissa of right side of the model
nx = 120 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
absorbright = .false.
absorbtop = .false.
absorbleft = .false.
# define the different regions of the model in the (nx,nz) spectral element mesh
nbregions = 10 # nb of regions and model number for each
1 120 1 120 1
1 120 121 132 2
1 120 133 144 3
1 120 145 156 4
1 120 157 168 5
1 120 169 180 6
1 120 181 192 7
1 120 193 204 8
1 120 205 216 9
1 120 217 222 10
-------------- next part --------------
#source 1. The components of a moment tensor source must be given in N.m, not in dyne.cm as in the DATA/CMTSOLUTION source file of the 3D version of the code.
source_surf = .false. # source inside the medium or at the surface
xs = 6500. # source location x in meters
zs = 5023. # source location z in meters
source_type = 2 # elastic force or acoustic pressure = 1 or moment tensor = 2
time_function_type = 3 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10.0 # dominant source frequency (Hz) if not Dirac or Heaviside
tshift = 0.0 # time shift when multi sources (if one source, must be zero)
anglesource = 0.0 # angle of the source (for a force only)
Mxx = -0.83 # Mxx component (for a moment tensor source only)
Mzz = -0.018 # Mzz component (for a moment tensor source only)
Mxz = -0.129 # Mxz component (for a moment tensor source only)
factor = 1.d0 # amplification factor
More information about the CIG-SEISMO
mailing list