forums.silverfrost.com Forum Index forums.silverfrost.com
Welcome to the Silverfrost forums
 
 FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

Numerical difference Silverfrost v8 vs Salford 2.54
Goto page Previous  1, 2
 
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Support
View previous topic :: View next topic  
Author Message
jcherw



Joined: 27 Sep 2018
Posts: 49
Location: Australia

PostPosted: Mon Mar 04, 2019 8:28 am    Post subject: Reply with quote

I finally managed to track down the issue after some further discussion with mecej4, who kindly recompiled the code.

It turns out the the version of Salford v2.54 which was passed on to me as a straight folder under c:\Program Files (x86) (ie not an install) compiled by default with double precision reals. I am not sure whether this is standard or influenced by the original setup. I am also not sure if the previous developed was aware of this. I had inherited BAT compile files from the previous developer which did not stipulate 'all double precision' compilation flags.

Silverfrost v8.3 compiles by default with single precision reals.

For most of the test problems I had at hand this did not cause any issue other than totally insignificant differences.

However, when I further analysed the test problem I posted (and that created the problem) did change the basic input via a unit conversion and single precision truncation significantly reduced the heat influx.

I am still not sure if the code require requires a blanket double precision for all other calcs, but obviously as a matter of safety it can't hurt. Most codes of this style I have worked with before had the double precision declarations explicitly inside the code, obviously to avoid the problem I ran into.

Lessons learned (for the umpteenth time): never make assumptions about defaults, and understanding the physics (which was pointed out to me by mecej4) is very valuable when tracking down bug.

All comments from various contributors to this thread are much appreciated.
Back to top
View user's profile Send private message
PaulLaidler
Site Admin


Joined: 21 Feb 2005
Posts: 5792
Location: Salford, UK

PostPosted: Mon Mar 04, 2019 8:57 am    Post subject: Reply with quote

These things are easy to miss. The original compilations may have used /config to configure the compiler to use /dreal by default.
Back to top
View user's profile Send private message
jcherw



Joined: 27 Sep 2018
Posts: 49
Location: Australia

PostPosted: Mon Mar 04, 2019 10:31 am    Post subject: Reply with quote

Is there a standard fortran statement that will declare all reals as kind=2 (double precision) ?
Back to top
View user's profile Send private message
PaulLaidler
Site Admin


Joined: 21 Feb 2005
Posts: 5792
Location: Salford, UK

PostPosted: Mon Mar 04, 2019 10:59 am    Post subject: Reply with quote

The safest way may be to call SELECTED_REAL_KIND to get the KIND value for the compiler.
Back to top
View user's profile Send private message
mecej4



Joined: 31 Oct 2006
Posts: 1084

PostPosted: Mon Mar 04, 2019 2:01 pm    Post subject: Re: Reply with quote

jcherw wrote:
Is there a standard fortran statement that will declare all reals as kind=2 (double precision) ?

First of all, note that kind numbers are NON-PORTABLE.

For real variables, the default kind numbers in FTN95 are 1 (32-bit real) and 2 (64-bit real). For FTN95, there is also the X87 10-byte real, for KIND=3. FTN95 provides the /ALT_KINDS option to enable the use of code that assumes kind numbers of 4 and 8 for 32-bit and 64-bit reals (with some limitation). Details may be seen in the FTN95 /help output.

As Paul wrote, in portable code you would use SELECTED_INT_KIND. For a project made up of a large number of source files, a somewhat shorter (but still portable) way is to put the following code into a separate module source file:
Code:
module floats
integer, parameter :: sp = kind(0.0),  dp = kind(0.0d0)
end module floats

In each subprogram of your project, you may add USE FLOATS and then change variable declarations from REAL to REAL(sp), if you choose single precision, or to REAL(dp) if you choose double precision.

If you choose this approach, where you have real constants in the code you may need to append "_sp" or "_dp". You have the value of Pi specified in various places in the code, but you do not have enough digits to provide double precision for Pi.
Back to top
View user's profile Send private message
jcherw



Joined: 27 Sep 2018
Posts: 49
Location: Australia

PostPosted: Tue Mar 05, 2019 1:26 am    Post subject: Reply with quote

Just for completeness, here is a small test program that clarifies the issue.

Quote:
Program testPrec

! Test precision of time conversion
! Heat flux (utu) unit coversion (/day to /sec) to comply with SI units
! which is used in remainder of program (solver etc)

implicit none

real*4 uqh4, utu4, res4 ! heat flux, conversion & result single precision
real*8 uqh8, utu8, res8 ! heat flux, conversion & result double precision
real*8 res48, dif8, tot8

utu4 = 86400.
utu8 = 86400.

uqh4 = -.1036E+09
uqh8 = -.1036E+09

res4 = uqh4/utu4
res8 = uqh8/utu8

write (*,'(A25,F15.5, E15.5, F15.5)') &
' Single utu4 uqh4 res4: ', utu4, uqh4, res4
write (*,'(A25,F15.5, E15.5, F15.5)') &
' Double utu8 uqh8 res8: ', utu8, uqh8, res8

res48 = res4
dif8 = res48-res8

write (*,'(A24,1X,3F15.5)') ' res48 res8 dif8:', res48, res8, dif8

! Accumulated heat flux difference after 1000 days (run of simulation)

tot8 = 1.D03*utu8*dif8
write (*,'(1X,A45,9X,F15.5)') &
' Accumulated heat flux difference @ 1000 days:', tot8

end

and the results

Quote:
Single utu4 uqh4 res4: 86400.00000 -0.10360E+09 -1199.07410
Double utu8 uqh8 res8: 86400.00000 -0.10360E+09 -1199.07407
res48 res8 dif8: -1199.07410 -1199.07407 -0.00002
Accumulated heat flux difference @ 1000 days -1953.12499
Back to top
View user's profile Send private message
mecej4



Joined: 31 Oct 2006
Posts: 1084

PostPosted: Tue Mar 05, 2019 3:30 am    Post subject: Reply with quote

The expected precision of 32-bit reals is seven significant decimal digits. You are taking the difference between two floating point numbers that differ only in their eighth digits; that difference in Q (between res48 and res8), as reported, is 2 or 3 X 10-6, in whatever physical units were used. You are then multiplying that by a large number (1000 days X 86400 seconds/day). The result of doing so, which is about 2,000, may appear large, but that is simply because of the units used. The fractional error between res48 and res8 is still about 2 X 10-8, and should not cause any worries.

Last edited by mecej4 on Thu Mar 07, 2019 4:48 am; edited 1 time in total
Back to top
View user's profile Send private message
LitusSaxonicum



Joined: 23 Aug 2005
Posts: 1960
Location: Yateley, Hants, UK

PostPosted: Wed Mar 06, 2019 5:09 pm    Post subject: Reply with quote

You don't need a completely standard way of doing anything once you locked yourself into FTNxx by using Clearwin+ - you just need a way that works, and tells you explicitly just in case SF changes the rules.

The best way of doing this is not through a compiler setting or command line option, but with an OPTIONS directive in every file.

Code:
      OPTIONS (DREAL)


If you try to compile this source code with a different compiler, it will baulk at the line.

You might even help someone in future by adding a comment:

Code:
      OPTIONS (DREAL) ! This makes all REAL variables DOUBLE PRECISION in FTN95 (non-standard)


Who knows, you might even add comments to further describe the issue solved:

Code:
      ! The right answers are not obtained with low precision REALs


For mecej4, I suspect that the round-off error was cumulative, as the results shown in the graph slowly diverge. It's that, rather than the 'big error' at the end of the procedure that count.

Eddie
Back to top
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Support All times are GMT + 1 Hour
Goto page Previous  1, 2
Page 2 of 2

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © 2001, 2005 phpBB Group