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

Mon May 7 06:56:41 PDT 2018

Hi Charles,

Perfect!
You are welcome.

Best regards,
Dimitri.

On 05/07/2018 03:55 PM, Chunyang Xiao wrote:
> Hi Dimitri,
> 
> 
> Thank you for your reply. I checked the line of 'Max stability for P 
> wave velocity...' as you told me to, and found it was over 0.5.
> 
> 
> After I lowered my time step, the problem was solved! 😊  Your 
> suggestions are really helpful! Thank you for your assistance.
> 
> 
> Best regards,
> 
> Charles
> 
> 
> ------------------------------------------------------------------------
> *From:* Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
> *Sent:* Friday, 4 May 2018 5:14:25 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 Charles,
> 
> Thanks! That is good news. Now it is probably simply that your time step
> is too high and you are above the CFL stability limit. Could you check
> what the following line says:
> 
> Max stability for P wave velocity (must be below about 0.50 or so) = ...
> 
> and check that the number is below 0.5?
> 
> Best regards,
> Dimitri.
> 
> On 05/03/2018 05:18 AM, Chunyang Xiao wrote:
>> 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
>> 
>> 
>> _______________________________________________
>> 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

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