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

[Chunyang Xiao]

Hi Dimitri,


Thank you so much for identifying the problem in my code. I didn't realize that the error message '0x4192BA in read_regions_ at read_regions.f90:99' means that I have problems in line 99 of read_regions.f90, and thank you for telling me.

I have exactly the error you said: I set Vp = 1, and Vs = 1 in the example Rayleigh_wave_no_crack. After I set Vp = 2, the error was gone.


Now as I changed my materials to be


1 1 2000.d0 3000.d0 2000.d0 0 0 9999.d0 9999.d0 0 0 0 0 0 0
2 1 1.d0 3000.d0 0.d0 0 0 9999.d0 9999.d0 0 0 0 0 0 0


An error came again, which said


#0  0x7F5BFED14E08
#1  0x7F5BFED13F90
#2  0x7F5BFE4454AF
#3  0x4489B6 in compute_forces_viscoelastic_ at compute_forces_viscoelastic.F90:601
#4  0x44E240 in compute_forces_viscoelastic_main_ at compute_forces_viscoelastic_calling_routine.F90:67
#5  0x4936C7 in iterate_time_ at iterate_time.F90:165
#6  0x40230A in specfem2d at specfem2D.F90:381
#7  0x7F5BFE43082F


Line 601 of compute_forces_viscoelastic.F90 says 'tempx1(i,j) = jacobianl * (sigma_xx*xixl + sigma_zx*xizl) ! this goes to accel_x'. I couldn't figure out what that means. I am thankful for any input or suggestion about it.


I am also curious about how to find max stability and output_solver.txt because floating point errors sometimes also happened when I changed the time step.

Thank you again for your assistance.

Sincerely,
Charles Xiao
________________________________
From: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Sent: Thursday, 3 May 2018 7:24 AM
To: cig-seismo at geodynamics.org; Chunyang Xiao
Subject: Re: [CIG-SEISMO] How to choose time step and other parameters to avoid floating-point error


Hi all,

I have added the stop statement.

Best regards,
Dimitri.

On 05/02/2018 11:20 PM, Dimitri Komatitsch wrote:
>
> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/cig-seismo/attachments/20180503/6050bc9b/attachment-0001.html>