Re: [fluka-discuss]: Numerical precision in FLUKA

From: Paola Sala <paola.sala_at_mi.infn.it>
Date: Sat, 17 Oct 2015 17:07:46 +0200

Dear Roberta
apparently you corrected only in one of the many statements used for
dumping the variables.
The lines that you sent come from the instruction
                    WRITE(*,200)'333333',XTUP(1),XTUP(2),
    $ XTUP(6),XTUP(7),XTUP(8)
where the xtups are REAL variables. Therefore the values ar rounded to
some "real precision.
I changed this to
WRITE(*,200)'333333',JTRACK,ETRACK-AM(JTRACK),
     $ CXTRCK,CYTRCK,CZTRCK
and I get the lines copied below.

Besides that:
-you cannot use the same "200" format for the two write statements, since
Jtrack is integer but xtup(1) is real. On g77 it might work, on gfortran
it crashes (correctly)
- in the IF clause on the forward angles, please use double precision
constants to compare with double precision variables
- Add a check on jtrack>0 to avoid scoring heavy fragments
(jtrack<0) for which am(jtrack) is not defined
- in the input, you do not need to define BEAMAXIS if the standard ones
(beam along z) are fne for you. If you wish to change, instead, check the
alignment of the card, as it is now the first what is read as 0 insttead
of 1.0

You do not need to OPEN the unit 42 in the input if you open it in mgdraw.

Hope ut helps
Paola


333333 1 0.398923488468048E+03
-0.787887208194129E-04 -0.350455774896683E-04 0.999999996282072E+00
 333333 1 0.398595058652988E+03
-0.144492666942190E-04 0.466919353024614E-04 0.999999998805541E+00
 333333 1 0.398904170469407E+03
0.308129771260977E-04 0.128930462558836E-05 0.999999999524449E+00
 333333 1 0.398893560709029E+03
-0.256446925081398E-04 0.573916701910749E-04 0.999999998024273E+00
 333333 1 0.398915938216952E+03
0.157569500918632E-03 0.106517493708166E-03 0.999999981912938E+00
 333333 1 0.398918559668654E+03
0.272263035060327E-04 -0.263461009859279E-04 0.999999999282306E+00
 333333 1 0.398895042319447E+03
0.491992483453783E-04 -0.176947858739383E-04 0.999999998633164E+00




> Dear Paola,
> I have corrected the statement, but I get the same issue.
> I am using lxplus, the fortran version is:
> GNU Fortran (GCC) 4.4.7 20120313 (Red Hat 4.4.7-16)
> and I am using FLUKA for 32bit.
> Attached are the input file and the mgdraw.f file.
> Thanks
> Roberta
>
>
>
> 2015-10-16 9:21 GMT+02:00 paola sala <paola.sala_at_cern.ch>:
>
>> Dear Roberta,
>> looking better at your message:
>> the write statement is not compatible with the format, Jtrack is an
>> integer, not a real. Also, the output unit (42 in your case) must be an
>> integer number. Depending on the compiler, these instructions can
>> trigger a
>> crash, or output a nonsense number for Jtrack, or... could you please
>> specify which compiler did you use, and, better correct the statements,
>> for
>> instance
>>
>> 200 FORMAT(1x,A,I5,4E25.15)
>>
>> WRITE(42,200)'333333',JTRACK,ETRACK-AM(JTRACK),
>> $ CXTRCK,CYTRCK,CZTRCK
>>
>>
>> Paola
>>
>> On 10/16/2015 07:23 AM, Paola Sala wrote:
>>
>>> Dear Roberta
>>> strange. Could you please send me the input and the user routines?
>>> Thanks
>>> Paola
>>>
>>>> Dear Paola,
>>>> thank for your suggestion, indeed in my modified mgdraw.f I was using
>>>> an
>>>> unfitting format in writing the txt file.
>>>>
>>>> Anyway, now, I fixed it, using:
>>>> 200 FORMAT(' ',A,5E25.15)
>>>> ......
>>>> WRITE(42.0,200),'333333',JTRACK,ETRACK-AM(JTRACK),
>>>> $ CXTRCK,CYTRCK,CZTRCK
>>>>
>>>> And the values I get for CZTRCK (for angle < 0.4mrad ) are:
>>>> CZTRCK=1 or CZTRCK=0.999999940395355E+00
>>>> and nothing between these values.
>>>>
>>>> See, as an example for angle grater than 0 and less than 0.4mrad, few
>>>> output lines below (*).
>>>>
>>>> So, now it doesn't seem a numerical precision issue anymore.
>>>> Do you have any idea on how to cope with that?
>>>> Thank you,
>>>> Roberta
>>>>
>>>> (*)
>>>> 333333 0.700000000000000E+01 0.489577674865723E+02
>>>> -0.379900448024273E-03 0.734752029529773E-04
>>>> 0.999999940395355E+00
>>>> 333333 0.100000000000000E+01 0.398913940429688E+03
>>>> 0.243720045546070E-03 0.187961297342554E-03
>>>> 0.999999940395355E+00
>>>> 333333 0.100000000000000E+01 0.378580841064453E+03
>>>> 0.240417313762009E-03 0.217997265281156E-03
>>>> 0.999999940395355E+00
>>>> 333333 0.100000000000000E+01 0.398872375488281E+03
>>>> -0.289755815174431E-03 -0.169058344908990E-03
>>>> 0.999999940395355E+00
>>>> 333333 0.100000000000000E+01 0.398769622802734E+03
>>>> -0.195283661014400E-03 -0.303217006148770E-03
>>>> 0.999999940395355E+00
>>>> 333333 0.100000000000000E+01 0.398898071289062E+03
>>>> 0.582441571168602E-04 0.355932890670374E-03
>>>> 0.999999940395355E+00
>>>> 333333 0.100000000000000E+01 0.398916778564453E+03
>>>> -0.868009956320748E-04 0.254166428931057E-03
>>>> 0.999999940395355E+00
>>>>
>>>> 2015-10-13 7:36 GMT+02:00 Paola Sala <paola.sala_at_mi.infn.it>:
>>>>
>>>> Hi
>>>>> it looks like there is some single precision variable, either in the
>>>>> writing or reading procedure, that cuts your rounding.
>>>>> Please check all the steps, and be sure to use double precision
>>>>> variables
>>>>> (the default in FLUKA) and to write them with enough decimals.
>>>>> You might also want to try the predefined FLUKA scoring options.
>>>>> Ciao
>>>>> Paola
>>>>>
>>>>>> Dear experts,
>>>>>> I have a question on the numerical precision which FLUKA can get on
>>>>>>
>>>>> the
>>>>>
>>>>>> angular distributions.
>>>>>> I simulated a 400 GeV proton beam incident on a thick Be target, 50
>>>>>> cm
>>>>>> long.
>>>>>> I modified the mgdraw.f routine in order to get an ASCII file with
>>>>>> the
>>>>>> information I need.
>>>>>> I am interested in the particles produced and going out of the
>>>>>> target
>>>>>>
>>>>> at
>>>>>
>>>>>> very small angle, say less than 0.2 mrad. I simulated my geometry
>>>>>> with
>>>>>> cylindrical symmetry along z axis, so I am interested in particles
>>>>>>
>>>>> with
>>>>>
>>>>>> ZSCO close to 1.
>>>>>> I see that the next value less than 1 is ZSCO=0.99999994,
>>>>>>
>>>>> corresponding
>>>>> to
>>>>>
>>>>>> an angle of 0.00034641 rad. So I cannot get a continuous
>>>>>> distribution
>>>>>>
>>>>> in
>>>>>
>>>>>> that angle range.
>>>>>> I am attaching two plots to show you the issue, they show the polar
>>>>>>
>>>>> angle
>>>>>
>>>>>> vs kinetic energy for neutrons, one is with angle less than 10 mrad,
>>>>>>
>>>>> and
>>>>> I
>>>>>
>>>>>> suppose that the distribution is reasonable, and the other is for
>>>>>>
>>>>> angle<1
>>>>>
>>>>>> mrad.
>>>>>> I am wondering if there is any way to increase the precision.
>>>>>> Also, can you please tell me if the distribution I get is the
>>>>>> expected
>>>>>>
>>>>> one
>>>>>
>>>>>> or if you see further issues with that?
>>>>>> Thank you,
>>>>>> best regards,
>>>>>> Roberta Volpe
>>>>>>
>>>>>>
>>>>> Paola Sala
>>>>> INFN Milano
>>>>> tel. Milano +39-0250317374
>>>>> tel. CERN +41-227679148
>>>>>
>>>>>
>>>>>
>>> Paola Sala
>>> INFN Milano
>>> tel. Milano +39-0250317374
>>> tel. CERN +41-227679148
>>>
>>> __________________________________________________________________________
>>> You can manage unsubscription from this mailing list at
>>> https://www.fluka.org/fluka.php?id=acc_info
>>>
>>>
>>
>


Paola Sala
INFN Milano
tel. Milano +39-0250317374
tel. CERN +41-227679148

__________________________________________________________________________
You can manage unsubscription from this mailing list at https://www.fluka.org/fluka.php?id=acc_info
Received on Sat Oct 17 2015 - 18:41:53 CEST

This archive was generated by hypermail 2.3.0 : Sat Oct 17 2015 - 18:41:54 CEST