CHALLENGE using root_fortran to find root of a complex function

Mubaraqlanre · September 8, 2023, 3:07pm

Greetings everyone

I am working on a modeling project that involves some complex computations. I transitioned a model from R to Fortran due to performance concerns that hindered the efficiency of the original model being to slow in R. Within this Fortran environment, I am in need of robust root-finding algorithms. While exploring, I came across the use of root-fortran library for finding root in modern Fortran. However on translating the r code to fortran and solving the code using similar tolerance as used in R code I have not been getting the same values as in the R code. The R function is a bit complex and translated to fortran as below

function Conc_calc(zHV) result(result)
real(Kind=8), intent(in) :: zHV
real(Kind=8) :: A, B, Kw, K1, K2, K3, H1, H2, zH
real(Kind=8) :: zH0, pH0, pHI, buff_L
real(kind=8), dimension(12) :: vec1, Vec2, vec3, vec4, vec5, vec6, vec7
real(kind=8), dimension(12) :: vec8, vec9, vec10, vec11, vec12
real(kind=8) :: Amatrix(12,12), Bvector(12)
real(kind=8) :: roots(12)
real(kind=8) :: result

        integer :: n=12 ! specify the size of the matrix
    
    
        A = x10 + x20 + x30 + InN - OutN
    
        B = y10 + y20 + y30 + y40 + InC
    
               
    
        Kw = 10.0**(-14.0)   ! Kw = 10.0**(-14.0) at 25 C
    
    
        !Print *, 'zHV', Tf, ti
                                    
        K1 = 5.67 * (10.0**(-10.0)) * exp(-6286.0 * ((1.0 / Tf) - (1.0 / 298.15))) !defines the first H2CO3 dissociation into HCO3-
    
        K2 = (10.0**(-1.0 * (((2902.39 / Tf) + (0.02379 * Tf) - 6.4980)))) !defines the second HCO3- dissociation into CO3--
    
        K3 = (10.0**(-1.0 * (((3404.71 / Tf) + (0.032786 * Tf) - 14.8435)))) !defines the third H2CO3-- dissociation into H2O and co2                   
    
        H1 = (56.0 * exp(4092.0 * ((1.0 / Tf) - (1.0 / 298.15)))) / ((1.0 / 8.314) * 0.001 * (1.0 / Tf) * &
    
                (1.013 * 10.0 ** 5)) 
    
        H2 = ((1.0 / 29.41) * exp(2400.0 * ((1.0 / Tf) - (1.0 / 298.15)))) / ((1.0 / 8.314) * 0.001 * &
    
                        (1.0 / Tf) * (1.013 * 10.0 ** 5))

                                      
    
        zH = zHV * V

        if(ti == 2 .and. ti3 == 3) then
            write(*,*) 'zH is ',zH, ti, ti3                        
        end if  
    
        vec1 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, -zH, 0.0_wp, & 
        0.0_wp, (1.0_wp + H1 * (V / Vg)), 0.0_wp/) ! x1
    
        vec2 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, K1* V, 0.0_wp, 0.0_wp, H1 * (V / Vg), 0.0_wp/) ! x2
        vec6 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
        vec3 = (/0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, -zH, 0.0_wp, 0.0_wp, H2* (V / Vg)/) ! y1
    
        vec4 = (/0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, K2*V, -zH, 0.0_wp, H2 * (V / Vg)/) ! y2
    
        vec5 = (/1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &  
        K3* V, 0.0_wp, (1.0_wp + H2 * (V / Vg))/) ! y3
    
        vec7 = (/0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! y4
    
        vec8 = (/1.0_wp, 0.0_wp, 0.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RH2
    
        vec9 = (/-1.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK3
    
        vec10 = (/0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK2
    
        vec11 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK1
    
        vec12 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RH1
    
                
    
        Amatrix = reshape([vec1, vec2, vec6, vec3, vec4, vec5, vec7, vec8, vec9, vec10, &
    
        vec11, vec12], shape(Amatrix))
    
                        ! Set the column names
    
        !character(2) :: colnames(12) = ["x1", "x2", "x3", "y1", "y2", "y3", "y4", "RH2", "RK3", &
    
                        ! "RK2", "RK1", "RH1"]
    
        Bvector = (/y30, y20 + InC, y10, y40 - OutC, x20 + InN, x10, x30 - OutN, 0.0_wp, 0.0_wp, 0.0_wp, &
        H1 * A * (V / Vg), H2 * B * (V / Vg)/) ! const
    
        call solve_system(n, Amatrix, Bvector, roots)   !I called lapack to solve the system         
    
        zH0 = zHV0 * V
        pH0 = -log10(zH0 / V)
        pHI = -log10(zH / V) !I renamed the line pH to PHI
        result =(zH0 - InC + (-roots(9) + roots(10) + roots(11)) + (buff_L * (pHI - pH0))) - ((zH))
    end function Conc_calc

I used the root module of root fortran and I called it as below

call root_scalar(‘brent’,Conc_calc, 0.0_wp, 1.0_wp, zHV, result, iflag, atol=1.0_wp,Ftol = 0.10_wp**(0.25_wp))

However for the last 3 days I have been stuck as I havent been getting the same result as the R code. I was wondering if I got a step wrong or not so I started to print every single variable. I noticed that when I print zH which is within the function I am finding the root of in both fortran and R I get the value to be in fortran as below

.\my_program.exe zH is0.0000000000000000
zH is 57.641325192598238
zH is 28.820662596299119
zH is 14.410331298149559
zH is 7.2051656490747797
zH is 3.6025828245373899
zH is 1.8012914122686949
zH is 0.90064570613434747
zH is 0.45032285306717373
zH is 0.22516142653358687
zH is 0.11258071326679343
zH is 5.6290356633396717E-002
zH is 2.8145178316698358E-002
zH is 1.4072589158349179E-002
zH is 7.0362945791745896E-003
zH is 3.5181472895872948E-003
zH is 1.7590736447936474E-003 the last line in fortran

where as in the R code I got values lower than this

1] zH is 0"
[1] " zH is 57.6413255382803"
[1] " zH is 28.8206627691402"
[1] " zH is 14.4103313845701"
[1] " zH is 7.20516569228504"
[1] " zH is 3.60258284614252"
[1] " zH is 1.80129142307126"
[1] " zH is 0.90064571153563"
[1] " zH is 0.450322855767815"
[1] " zH is 0.225161427883908"
[1] " zH is 0.112580713941954"
[1] " zH is 0.0562903569709769"
[1] " zH is 0.0281451784854884"
[1] " zH is 0.0140725892427442"
[1] " zH is 0.00703629462137211"
[1] " zH is 0.00351814731068605"
[1] " zH is 0.00175907365534303" where the fortran ended
[1] " zH is 0.000879536827671514"
[1] " zH is 0.000489560168170512"
[1] " zH is 0.000674886029990813"
[1] " zH is 0.000650298707008735"
[1] " zH is 0.00064673969984793"
[1] " zH is 0.000646746327972608"
[1] " zH is 0.000646746309843868"
[1] " zH is 0.000646746309843775"
[1] " zH is 0.000646746309843776"
[1] " zH is 0.000646746309843776"

After using the one specified in the R code since in R, the default setting is. tol = .Machine$double.eps^0.25.
I have really been tweaking the ftol value which has been changing the values but stopped at the value above in the call rootscalar**(used 1.0_wp**(0.25_wp)** )in fortran. I am working with ecological data and this step is an important one to finish up with the translation.

I would be very grateful if I get suggestions on the next step to take to figure out this as I have been stuck on this step for almost a week. Many thanks in advance

FortranFan · September 8, 2023, 3:29pm

@Mubaraqlanre ,

Sorry to read of your trouble, hopefully you will able to resolve it quickly now.

I barely spent a few seconds at your problem and the first glace was at an equation and a thought immediately came to mind you need to apply consistent precision between R and Fortran in your arithmetic expressions involving literal constants if you want consistent results between the two.

Chances are R being effectively an interpreted environment based on GCC ecosystem uses the C language default of double for floating-point literal constants.
Whereas Fortran given its legacy dating back to mid-1950s uses default real which for your compiler might be the equivalent of C float.

I suggest you check the precision in R and if it is, as I suspect, the equivalent of C double, then reauthor your Fortran code suitably with the use of:

defined KIND e.g., use, intrinsic :: iso_c_binding, only : WP => c_double,
use WP as the defined KIND consistenly in all floating-point object and constant declarations in Fortran - see point #3 below with TREF.
use named constants for enhanced mnemonics in your formulae e.g., real(WP), parameter :: TREF = 298.15_wp to define the literal constants only once,
use the constants such as TREF and INVTREF in the Fortran equations.

Mubaraqlanre · September 8, 2023, 4:14pm

Many thanks for the suggestions and inputs, really appreciate

For 1. I did define the kind and the function is within a subroutine where Tf changing so not a constant. I used the example provided on root fortran sample online. is there need to change this?

subroutine pH_calculation
use root_module, only: root_scalar, wp => root_module_rk

Regarding the 2. Now I have adjusted the wp instead of kind=8 but I still get the same values

For 3 and 4, I do have parameters defined as the code is part of a bigger model and virtually all values are the same with the R code except for this aspect for the subroutine. For Tf in that line was used to refer to a variable within the subroutine that is related to the function, so was not a constant.
The subroutine is very long but I will share part of it below

subroutine pH_calculation
use root_module, only: root_scalar, wp => root_module_rk
character(len=5), parameter :: method = ‘brent’
integer :: ti3, ti
real(wp) :: zHV, Tf, V, Vg, x10, x20, x30, y10, y20, y30
real (wp) :: y40, InC, InCf, InN, OutN, zHV0, OutC
real(wp) :: outnf, InNf, outcf, Zh0f
real(wp) :: Tfix, Volumefix, volumegfix, x10f
real(wp) :: x20f, x30f, y10f, y20f, y30f, y40f
integer:: iflag
!real(kind=8), dimension(12) :: roots
real(wp) :: result
real(wp) :: roots(12)

        Do ti3=1, cv 
            Do ti = 2,tn
                if (DEBUG_pH) then ! pH module
    
                    if (ti == 2) then
    
                    y10f = 0
    
                    y20f = 0
    
                    y30f = 0
    
                    y40f = 0
    
                    x10f = 0
    
                    x20f = 0
    
                    x30f = 0
                end if
    
                InCf = TANproduction(ti, ti3) * (ASL/WSL) / 14.0 / 2.0
                
                OutNf = (SoilEmissionFlux(ti-1, ti3) * (PatchArea_L/100.0)) / 14.0 + (((TANBudget(ti-1, ti3)) * &
                        (1.0-LOSS)) / 14.0)

                InNf = TANproduction(ti, ti3) * (ASL/WSL) / 14.0

                Tfix = Ts(ti) + 273.15 ! K

                Volumefix = WaterBudget2cm(ti, ti3) * (ASL/WSL) ! dm3
                Volumegfix = (Porosity * PatchArea_L * ASL) - Volumefix ! dm3

                OutCf = 0.0
               
                zH0f = 10.0 ** (-pH(ti-1, ti3)) !Issue with this line
                

                Tf=Tfix !To be passed to the function
                
                V=Volumefix
            
                Vg=Volumegfix
            
                x10=x10f 
            
                x20=x20f
            
                x30=x30f
            
                y10=y10f
            
                y20=y20f
            
                y30=y30f
            
                y40=y40f
            
                InC=InCf
            
                InN=InNf
            
                OutN=OutNf
            
                zHV0=zH0f
            
                OutC=OutCf
                !print *, 'Vg is',  InCf,  ti, ti3
                !print *, 'outN', zH0f, ti, ti3
                call root_scalar('brent',Conc_calc, 0.0_wp, 1.0_wp, zHV, result, iflag, atol=1.0_wp,Ftol = 0.10_wp**(0.25_wp))
                
                Hconc(ti, ti3)= zHV ! the root of the function is the H+ concentration
                !print *, 'Vg', ZHV, ti
                pH(ti, ti3)= -log10(Hconc(ti, ti3))
                 
                roots =Conc_calc2(zHV)
                y10f = roots(4)
                y20f = roots(5)
                y30f = roots(6)
                y40f = roots(7)
                x10f = roots(1)
                x20f = roots(2) 
                x30f =roots(3)
                
                !test = reshape([roots(8:12), InC], shape(test))
              
            else
                pH(ti, ti3) = constpH
                Hconc(ti, ti3) = 10.0 ** (-pH(ti,ti3))
            end if
       
        end do 
    end do 
    contains 
    function Conc_calc(zHV) result(result)
        real(wp), intent(in)  :: zHV
        real(wp) :: A, B, Kw, K1, K2, K3, H1, H2, zH
        real(wp) :: zH0, pH0, pHI, buff_L
        real(wp), dimension(12) :: vec1, Vec2, vec3, vec4, vec5, vec6, vec7
        real(wp), dimension(12) :: vec8, vec9, vec10, vec11, vec12
        real(wp) :: Amatrix(12,12), Bvector(12) 
        real(wp) :: roots(12)
        real(wp) :: result
        
        integer :: n=12 ! specify the size of the matrix
    
    
        A = x10 + x20 + x30 + InN - OutN
    
        B = y10 + y20 + y30 + y40 + InC
    
               
    
        Kw = 10.0**(-14.0)   ! Kw = 10.0**(-14.0) at 25 C
    
    
        !Print *, 'zHV', Tf, ti
                                    
        K1 = 5.67 * (10.0**(-10.0)) * exp(-6286.0 * ((1.0 / Tf) - (1.0 / 298.15))) !defines the first H2CO3 dissociation into HCO3-
    
        K2 = (10.0**(-1.0 * (((2902.39 / Tf) + (0.02379 * Tf) - 6.4980)))) !defines the second HCO3- dissociation into CO3--
    
        K3 = (10.0**(-1.0 * (((3404.71 / Tf) + (0.032786 * Tf) - 14.8435)))) !defines the third H2CO3-- dissociation into H2O and co2                   
    
        H1 = (56.0 * exp(4092.0 * ((1.0 / Tf) - (1.0 / 298.15)))) / ((1.0 / 8.314) * 0.001 * (1.0 / Tf) * &
    
                (1.013 * 10.0 ** 5)) 
    
        H2 = ((1.0 / 29.41) * exp(2400.0 * ((1.0 / Tf) - (1.0 / 298.15)))) / ((1.0 / 8.314) * 0.001 * &
    
                        (1.0 / Tf) * (1.013 * 10.0 ** 5))

                                      
    
        zH = zHV * V
    
        vec1 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, -zH, 0.0_wp, & 
        0.0_wp, (1.0_wp + H1 * (V / Vg)), 0.0_wp/) ! x1
        
        if(ti == 2 .and. ti3 == 3) then
            write(*,*) 'vec1 is ', zH             
        end if 
        
        vec2 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, K1* V, 0.0_wp, 0.0_wp, H1 * (V / Vg), 0.0_wp/) ! x2
        vec6 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
        vec3 = (/0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, -zH, 0.0_wp, 0.0_wp, H2* (V / Vg)/) ! y1
    
        vec4 = (/0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, K2*V, -zH, 0.0_wp, H2 * (V / Vg)/) ! y2
    
        vec5 = (/1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &  
        K3* V, 0.0_wp, (1.0_wp + H2 * (V / Vg))/) ! y3
    
        vec7 = (/0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! y4
    
        vec8 = (/1.0_wp, 0.0_wp, 0.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RH2
    
        vec9 = (/-1.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK3
    
        vec10 = (/0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK2
    
        vec11 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK1
    
        vec12 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RH1
    
                
    
        Amatrix = reshape([vec1, vec2, vec6, vec3, vec4, vec5, vec7, vec8, vec9, vec10, &
    
        vec11, vec12], shape(Amatrix))
    
                        ! Set the column names
    
        !character(2) :: colnames(12) = ["x1", "x2", "x3", "y1", "y2", "y3", "y4", "RH2", "RK3", &
    
                        ! "RK2", "RK1", "RH1"]
    
        Bvector = (/y30, y20 + InC, y10, y40 - OutC, x20 + InN, x10, x30 - OutN, 0.0_wp, 0.0_wp, 0.0_wp, &
        H1 * A * (V / Vg), H2 * B * (V / Vg)/) ! const
    
        call solve_system(n, Amatrix, Bvector, roots)   !I called lapack to solve the system         
    
        zH0 = zHV0 * V
        pH0 = -log10(zH0 / V)
        pHI = -log10(zH / V) !I renamed the line pH to PHI
        result =(zH0 - InC + (-roots(9) + roots(10) + roots(11)) + (buff_L * (pHI - pH0))) - ((zH))
    end function Conc_calc

end subroutine

Thanks a lot for the feedback but could be issue with ftol in root fortran?

Beliavsky · September 8, 2023, 5:48pm

A recent arXiv preprint Finding roots of complex analytic functions via generalized colleague matrices by Zhang and Rokhlin says

We implemented our algorithms in FORTRAN 77, and compiled it using GNU Fortran,
version 12.2.0.

I don’t see the code online, but the authors may be willing to send it if asked.

Mubaraqlanre · September 8, 2023, 7:27pm

Thanks a lot for the recommendation.
I will do a mail to the authors as soon as possible

FortranFan · September 8, 2023, 8:22pm

@Mubaraqlanre ,

I suggest first setting up your formulae securely in your Fortran code in a way that will be also consistent with the calculations in R.

Consider your equation toward what appears to be enthalpy of second reaction (or is it with Henry’s law constant) that has a lead constant. Your Fortran code has it as 1.0 / 29.41. You will realize that will result in a difference in the computations in R, likely along the lines shown below:

   blk1: block
      real(kind=8) :: a
      a = 1.0 / 29.41
      print *, "Block 1: default precision"
      print "(*(g0))", "a = 1.0/29.41: ", a
      print "(g0,b0)", "a (binary) : ", a
   end block blk1
   print *
   blk2: block
      use, intrinsic :: iso_c_binding, only : WP => c_double
      real(kind=WP), parameter :: ONE = 1.0_wp
      real(kind=WP), parameter :: x = 29.41_wp
      real(kind=WP), parameter :: a = ONE / x
      print *, "Block 2: recommended approach"
      print "(*(g0))", "a = 1.0/29.41: ", a
      print "(g0,b0)", "a (binary) : ", a
   end block blk2
end

C:\temp>gfortran -ffree-form p.f -o p.exe

C:\temp>p.exe
 Block 1: default precision
a = 1.0/29.41: 0.34002039581537247E-1
a (binary) : 11111110100001011010001011011100100000000000000000000000000000

 Block 2: recommended approach
a = 1.0/29.41: 0.34002040122407345E-1
a (binary) : 11111110100001011010001011011100100100101001010110001011001100

Note here block 1 is equivalent your code in Fortran whereas block 2 is likely equivalent to what occurs in R. Note the difference in the lead constant in the two blocks.

You are dealing with quantities with small numerical magnitudes as pertinent to aqueous ionic chemistry due to limited solubility of CO2 in water leading toward low concentrations of carbonic acid and bicarbonate in the aqueous phase.

So you want to first get the numerical precision consistent in your calculations and then you can proceed toward the solution of the equations and cross-checking the roots you get in Fortran against your R program (the latter albeit at a slower compute speed due to interpreted code at run-time).

And a way to get the numerical precision right in your Fortran code is to get the constants in your equations computed during compilation time with consistently defined KINDs as I show in block 2 above rather than simply writing them as 1.0 / 29.41 in Fortran.

Caveat emptor with floating-point literal constants in Fortran.

Harper · September 9, 2023, 12:04am

I saw Kw = 10.0**(-14.0) in your program. Why force the machine to do an unnecessary LOG and EXP? If 6-figure accuracy had been enough I would have written Kw = 1.0e-14 but as you were after higher accuracy I would use Kw = 1.0e-14_wp; others have told you about wp.

gareth · September 9, 2023, 4:23am

As others have stated above, in Fortran real constants are assumed single precision unless you give them a kind suffix – even if the LHS variables are double precision.

So in your code, a line like:

K1 = 5.67 * (10.0**(-10.0)) * exp(-6286.0 * ((1.0 / Tf) - (1.0 / 298.15)))

involves lots of single precision calculation on the RHS. So K1 will have less accuracy than you’d expect from a double precision calculation.

The solution is to provide a suffix every time you use a real constant, e.g.

K1 = 5.67_wp * (10.0_wp**(-10.0_wp)) * exp(-6286.0_wp * ((1.0_wp / Tf) - (1.0_wp / 298.15_wp)))

In Fortran you should always provide suffixes for real constants, to avoid such problems.

RonShepard · September 9, 2023, 7:25am

As pointed out elsewhere, the floating point exponents with integer values should almost always be written as integers. However, in this case, there is also the issue of writing the floating point literals in the usual way.

K1 = 5.67e-10_wp * exp(-6286.0_wp * ((1.0_wp / Tf) - (1.0_wp / 298.15_wp)))

Of course, it isn’t necessary to do that everywhere because fortran will probably do some of them correctly anyway because of the promotion rules for mixed-kind expressions, but I think it is a good idea for the programmer, especially a new fortran programmer, to specify the kinds explicitly anyway to make his intentions clear.

jacobwilliams · September 10, 2023, 3:16pm

Single precision strikes again!

It’s a rite of passage for every Fortran programmer to get bitten by this issue.

TL;DR

Add _wp to every floating point value, no matter how much you don’t want to. The syntax is 3.0_wp, 10.1e3_wp, etc.
Don’t use 3_wp expecting it to be a floating point value (if it doesn’t have the . then it is an integer).
Don’t use 3.0 without the _wp ever under any circumstances.
Don’t declare variables as real or real(8) or real*8 under any circumstances, always real(wp).
Define wp, for example using use iso_fortran_env, only: wp => real64 to get double precision. Put it in a module if you want to use it in multiple files.

FortranFan · September 10, 2023, 11:51pm

@Mubaraqlanre ,

Re: “transitioned a model from R to Fortran due to performance concerns … original model being too slow in R …,” please know you can speed up calculations significantly using Fortran (or, for that matter, C++ or C), basically any compiler based language that processes and optimizes your code directly to your compute platform. Fortran is really good at this once you get over a few of its legacy aspects.

Since your calculation involved the pH of the aqueous system of interest containing dissolved CO2, please revisit the system equations (and if you are up to it and can do so, to share them with the readers of this Discourse and any other programming / numerical solution / chemistry forums you follow) and consider:

formulating in terms of logarithms (natural log is, well, natural) that extends naturally with ln K the natural log of equilibrium constants and the expressions toward them in chemistry literature. You can then solve for the log of hydronium ion concentration(s). This may help you avoid most coding instructions toward exponentiation and the use of exp funciton and this can help you further with performance.
simple base transformation and/or arithmetic will then yield the pH for you.

FortranFan · September 11, 2023, 1:04am

So consider the K equation you have that “defines the first H2CO3 dissociation into HCO3-” which is in the form of the classic Arrhenius expression:

k = a \ exp(\frac { - E } { RT })

that is “scaled” in a sense against a reference temperature as

k = a \ exp(\ - E` [\frac { 1 } {T} - \frac { 1 } {T_{ref}}]\ )

which you seek to write in Fortran with literal constants as

K1 = 5.67 × 10^{-10} \ exp(\ - 6286.0 \ [\frac { 1 } {T} - \frac { 1 } {298.15}]\ )

So given where things stand with modern Fortran, please see the snippet below that contrasts what you show in the original post with an alternate approach consistent with what R likely calculates for you:

Click here for code

   blk1: block
      real(kind=8) :: K1
      real(kind=8) :: Tf
      print *, "Original post:"
      Tf = 400.0  !
      !defines the first H2CO3 dissociation into HCO3-
      K1 = 5.67 * (10.0**(-10.0)) * exp(-6286.0 * ((1.0 / Tf) - (1.0 / 298.15)))
      print "(*(g0))", "K1(decimal) for Tf = 400 K: ", K1
      print "(g0,b0)", "K1(decimal) for Tf = 400 K: ", K1
   end block blk1
   print *
   blk2: block
      use, intrinsic :: iso_c_binding, only : WP => c_double
      real(kind=WP) :: K1
      real(kind=WP) :: Tf
      real(kind=WP), parameter :: ONE = 1.0_wp
      real(kind=WP), parameter :: TREF = 298.15_wp
      real(kind=WP), parameter :: A1 = 5.67E-10_wp
      real(kind=WP), parameter :: E1 = 6286.0_wp
      print *, "Another approach, more secure in Fortran:"
      Tf = 400.0_wp  !
      !defines the first H2CO3 dissociation into HCO3-
      K1 = A1 * exp(-E1 * ( ONE/Tf - ONE/TREF ) )
      print "(*(g0))", "K1(decimal) for Tf = 400 K: ", K1
      print "(g0,b0)", "K1(decimal) for Tf = 400 K: ", K1
   end block blk2
end

Click here for compiler response

C:\temp>gfortran -ffree-form p.f -o p.exe

C:\temp>p.exe
 Original post:
K1(decimal) for Tf = 400 K: 0.12162608824641047E-6
K1(decimal) for Tf = 400 K: 11111010000000010100110000101001100010110110100011111110010100

 Another approach, more secure in Fortran:
K1(decimal) for Tf = 400 K: 0.12162606940969464E-6
K1(decimal) for Tf = 400 K: 11111010000000010100110000101000111000011011111001111011111000

Note you are dealing with a weak reaction and the absolute value of the equilibrium constant is small:

first, notice the difference in the result in the two approaches,
second, better to work with the natural log of such an equilibrium constant rather than K1 directly in your equations to be solved.

FortranFan · September 11, 2023, 1:26am

So better yet will be for you to write the scaled expression as

ln\ k = ln\ a\ - (\frac { E` } {T_{ref}})\ [ t_{r}\ - 1 ]

where

t_{r}\ = \frac { T_{ref} } { T }

which you can author in Fortran as follows and which can help boost your code performance:

Click here for code

      use, intrinsic :: iso_c_binding, only : WP => c_double
      real(kind=WP), parameter :: ONE = 1.0_wp
      real(kind=WP), parameter :: a1 = 5.67E-10_wp  ! units?
      real(kind=WP), parameter :: lna1 = log(a1)
      real(kind=WP), parameter :: E1 = 6286.0_wp    ! Activation energy, kcal/mol?
      real(kind=WP), parameter :: TREF = 298.15_wp  ! degrees Kelvin
      real(kind=WP), parameter :: E1TREF = E1 / TREF
      real(kind=WP) :: K1
      real(kind=WP) :: lnK1
      real(kind=WP) :: Tf
      real(kind=WP) :: Tr
      Tf = 400.0_wp
      Tr = TREF / Tf
      lnK1 = lna1 - E1TREF*( Tr - ONE )

Mubaraqlanre · September 11, 2023, 7:23am

Thanks to everyone. I have learnt a lot from all the valuable comments and have been adjusting my code appropriately. I have seen some improvements in the value I obtain. I will continue to make changes implementing all the suggestions. Many thanks to everyone for helping with advice.

pcosta · September 11, 2023, 10:46am

I thought that we won’t face problems when using numbers that can be perfectly described with single precision format, like 3.0…

! cat test.f90
print*,3.d0/acos(-1.d0),3.0/acos(-1.d0)
end

Running this program with a recent version of gfortran produces the following output:

0.95492965855137202       0.95492965855137202

So in cases like this in mixed precision mode, it is fine, right?

jacobwilliams · September 11, 2023, 1:13pm

Why risk it? Always use _wp and you will never have a problem. Also:

don’t use .d0 under any circumstances.

pcosta · September 11, 2023, 1:23pm

Thanks! Though I don’t think I’m risking anything here. I agree with avoiding .d0, was just for brevity of the code snippet.

ofmla · September 11, 2023, 1:23pm

within the context of this thread shouldn’t it be lna1 = dlog(a1) instead of lna1 =log(a1) ? in this case it seems that it doesn’t make difference but I guess that in other cases it does

ivanpribec · September 11, 2023, 1:47pm

I guess this is true for whole numbers, and in case your processor uses base-2, the set of dyadic rationals, i.e. the numbers which can be represented exactly as a fraction of the from i/2^j where i and j are integers. For example 1.5 is 3/2^1 as a fraction.

But on a processor with a different word-size, and different floating point format, the constant 1.5 might not be promoted correctly, unless you add the kind specifier suffix.

RonShepard · September 11, 2023, 4:46pm

Yes, and an experienced fortran programmer who understands how the fortran KIND system works, how floating point values of various KINDs are stored, and how mixed-type expressions are converted by default might write code like that. But this discussion is about a new fortran programmer just learning the language, so it is probably better to write out all constants and even all conversions explicitly rather than relying on default conversion conventions.

BTW, there is no difference in the compiled code. Whether you write x*n or x*real(n,kind(x)), the compiler must do exactly the same thing in both cases, and exactly the same optimizations are allowed. It isn’t like the first expression allows some kind of shortcut that the second expression forbids. And fortran is not an interpreted language, so the number of characters does not affect performance.

Topic		Replies	Views
Announcing Polyroots-Fortran Announcements	3	947	October 6, 2022
General roots of cubic polynomial with perturbed R	9	1139	November 20, 2020
How do I use the gsl_root_fsolver_bisection from gsl in Fortran Help	6	969	September 17, 2022
Optimization, root finding, and nonlinear equation solving Announcements	11	3635	July 29, 2022
Integer root of floating point value	25	1218	February 24, 2023

CHALLENGE using root_fortran to find root of a complex function

Related topics