[CIG-SEISMO] Time step value and CPML, displacement blowing up

Wed Mar 22 07:54:47 PDT 2017

Dear SPECFEM3D developers,

my name is Mathias Pilch and I’m a physics student in Germany and working on my 
master's thesis. I’m using SPECFEM3D to simulate seismic signals of tracked vehicles,
which is done for a verification project of the Physics and Disarmament Group of Lehrstuhl
Experimentelle Physik III at the Technische Universität Dortmund. I am sending
a copy of this mail to my advisor Juergen Altmann.

We find SPECFEM3D an extremely useful tool and are very grateful to you for
developing and maintaining it.

I’m currently using the devel branch of SPECFEM3D.

There seems to be a problem concerning the usage of CPML absorbing boundary
conditions and using a time step value DT that is below a certain value.
In my case, under many different conditions, there seems to be a threshold value of
1.75d-5 for DT when using CPML. Any value that is smaller than that leads to a crash of the 
output solver.

The models I’m using are a variation of 			

	EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_elastic_absorbing_CPML_5sides

with the difference being a smaller geometry, e. g. 120 m * 9 m * 12 m, with element size 0.2 m 
resulting in 600 * 45 * 60 = 1620000 elements. Basically I’m making sure that all element edges have
the same length.  

In this example, the specfem3D solver fails with 

	*** Minimum period resolved = 7.22559926E-04 
        *** Maximum suggested time step = 1.56957449E-05 

	*** for DT : 1.0000000000000001E-005 
 	*** Max stability for wave velocities = 0.318557680

even though the DT chosen is below dt_suggested, and cmax = max stability for wave velocities is
below 0.5. It also fails when taking a value closer to dt_suggested, e.g. 1.5d-5. See attachment 
"model_with_120*9*12" for output and log files. 

The usual error message when using a value below 1.75d-5 looks like:

	Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

	Backtrace for this error:

	Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

	Backtrace for this error:
	#0 0x2AF167374B97
	#1 0x2AF167373D90 
	#2 0x3D3D83269F 
	#0 0x2B3DD69DCB97 
	#1 0x2B3DD69DBD90 
	#2 0x3D3D83269F 
	#3 0x404DD1 in check_stability_ at check_stability.f90:91 (discriminator 13)
	#3 0x404DD1 in check_stability_ at check_stability.f90:91 (discriminator 13) 
	#4 0x478376 in iterate_time_ at iterate_time.F90:185 
	#4 0x478376 in iterate_time_ at iterate_time.F90:185 
	#5 0x4FA00D in specfem3d_ at specfem3D.f90:375 
	#6 0x40319F in xspecfem3d at program_specfem3D.f90:34 
	#7 0x3D3D81ED5C 
	#5 0x4FA00D in specfem3d_ at specfem3D.f90:375 
	#6 0x40319F in xspecfem3d at program_specfem3D.f90:34 
	#7 0x3D3D81ED5C
	[lido-bl43112][[46853,1],104][btl_tcp_frag.c:237:mca_btl_tcp_frag_recv] mca_btl_tcp_frag_recv: readv failed: Connection
	reset by peer (104)
	[lido-bl43141][[46853,1],101][btl_tcp_frag.c:237:mca_btl_tcp_frag_recv] mca_btl_tcp_frag_recv: readv failed: Connection
	reset by peer (104)
	-------------------------------------------------------------------------- 
	mpirun noticed that process rank 102 with PID 0 on node lido-bl43131.lidocluster.hp exited on signal 8 (Floating point
	exception).

Which is I think caused by the 

	Max norm displacement vector U in all slices (m)

blowing up when running the solver. 

I used the .f90 programs in utils/CPML/, namely

	xadd_CPML_layers_to_an_existing_mesh
	xconvert_external_layers_of_a_given_mesh_to_CPML_layers

to create the CPML layers with a thickness of three element sizes, which seems to work as expected.
Surprisingly the code will not crash when using a thickness of only one element size, in any other
case it crashes again.  I also varied the percentage of CPML elements in the mesh, whether 
it’s around 50 % or 20 %, it still crashes. I don’t think there is an upper bound on CPML percentage,
but I wanted to verify that.

To rule out any error on my part I took the model

	EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_elastic_absorbing_CPML_5sides

and set 

	DT = 1.d-5
	NSTEPS = 1000000.

in DATA/Par_file. Lo and behold, it crashes again. Please take a look at 

	max_norm_displacement.pdf

in the attachment directory “cpml_example_edited".

which shows quite a weird behaviour. The half duration (or f0) is set to 5.0,
so at -10 s there should be no signal, but still, the displacement blows up even at that point. 

Those crashes only occur when using CPML as boundary condition and when DT is too low,
however with Stacey it works fine with any DT value, although the absorption isn’t nearly as desirable.

Switching to the SPECFEM3D master branch doesn’t help, the same error persists. 

May it be possible that this a fundamental problem of CPML, when DT values are below 1.75d-5?
Could it be that another or a second stability condition applies for CPML?
Or could this still be a problem on my side, and it works completely as intended on yours?

Best regards and thank you for developing SPECFEM3D.
Mathias

-------------- next part --------------
A non-text attachment was scrubbed...
Name: model_with_120*9*12.zip
Type: application/zip
Size: 2312343 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/cig-seismo/attachments/20170322/db54e27b/attachment-0002.zip>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: cpml_example_edited.zip
Type: application/zip
Size: 1339463 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/cig-seismo/attachments/20170322/db54e27b/attachment-0003.zip>
-------------- next part --------------