[CIG-SEISMO] How to choose time step and other parameters to avoid floating-point error

[Dimitri Komatitsch]

Hi Charles,

Thanks for your message!
Line 99 of src/meshfem2D/read_regions.f90 is this:

        poisson_ratio = 0.5d0*(vpregion*vpregion - 
2.d0*vsregion*vsregion) / (vpregion*vpregion - vsregion*vsregion)

thus it is just an error in your input file (Par_file), you have at 
least one material for which Vs = Vp and thus you get a division by zero.

Materials always have a value of Poisson's ratio > -1 (in practice, for 
natural materials, >= 0), thus for all natural materials you must have 
Vs < Vp / sqrt(2).

I will add a stop statement in the code, right before that line, to 
detect that.

Best regards,
Dimitri.

On 05/01/2018 02:33 PM, Chunyang Xiao wrote:
> Hi there,
> 
> 
> I am an undergraduate student majoring in mechanical engineering. 
> Recently when I changed the speed of P-wave (Vp from the section 
> 'velocity and density models' in Par_file) in SPECFEM2D, an error occurred.
> 
> 
>     Program received signal SIGFPE: Floating-point exception - erroneous
>     arithmetic operation.
> 
>     Backtrace for this error:
> 
>     #0  0x7F209C4FDE08
> 
>     #1  0x7F209C4FCF90
> 
>     #2  0x7F209C14D4AF
> 
>     #3  0x4192BA in read_regions_ at read_regions.f90:99
> 
>     #4  0x418CEE in read_parameter_file_ at read_parameter_file.F90:88
> 
>     #5  0x42521C in MAIN__ at meshfem2D.F90:403
> 
>     Floating point exception (core dumped)
> 
> 
> This happened a few times also when changed other parameters, such as 
> time step.  A similar error happened to Mr. Moritz Wehler before, and 
> Dimitri's answer on Wed May 27 07:40:46 PDT 2015 was
> 
> 
>     It seems you simply forgot to also divide the time step by two to
>     remain
> 
>     stable when you increased mesh resolution by two:
> 
>        *** Max stability for P wave velocity (must be below about 0.50
>     or so)
> 
>     = 0.972936358376379
> 
> 
> I was curious where we could see the 'Max stability for P wave velocity' 
> for a given case. As I turned to the SPECFEM2D user manual (version 7.0) 
> for help, I saw chapter 4.6 (How to choose the time step) say:
> 
> 
>     You can directly compare these values with the value given in
>     sentence ‘_Max stability for P wave velocity_’ in file
>     _output_solver.txt _to see whether you set the correct ∆t in
>     Par_file or not. For elastic simulation, the CFL value given in
>     _output_solver.txt _does not consider the V p /V s ratio......
> 
> 
> That makes me more confused because I can find neither the sentence ‘Max 
> stability for P wave velocity’ nor the file output_solver.txt. The file 
> named output_solver.txt only exists in the folders of a few examples I 
> did not use or run.
> 
> Also, in Ms. Parvathy Sindhu Nair's question on Mon Feb 15 14:25:50 PST 
> 2016, she and Dimitri mentioned CFL, which seems to matter but does not 
> appear in my Par_file.
> 
>     this is either the time step that is a bit too high (unlikely, since
>     the CFL you used (0.44) seems fine), or if the amplitude of the
>     source you used is very high......
> 
> 
> Would anyone kindly let me know how I can choose appropriate parameters 
> (time step, speed, etc.) to avoid such floating-point errors and where I 
> can see my CFL, output_solver.txt, and max stability for P wave?
> 
> Sincerely,
> Charles Xiao
> 
> 
> _______________________________________________
> CIG-SEISMO mailing list
> CIG-SEISMO at geodynamics.org
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-seismo
> 

-- 
Dimitri Komatitsch, CNRS Research Director (DR CNRS)
Laboratory of Mechanics and Acoustics, Marseille, France
http://komatitsch.free.fr