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

From: roberta volpe <roberta.volpe_at_cern.ch>
Date: Sat, 17 Oct 2015 18:03:37 +0200

Dear Paola,
thanks a lot.
Indeed I forgot to fix the line I used for the check.
Now I managed to get very small angle values.
Thanks also for the other suggestions, I need also to change the beam axis,
but in that case I was using the right alignment.
Regards,
Roberta

2015-10-17 17:07 GMT+02:00 Paola Sala <paola.sala_at_mi.infn.it>:

> 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 - 19:21:26 CEST

This archive was generated by hypermail 2.3.0 : Sat Oct 17 2015 - 19:21:27 CEST