Thanks a lot for all your valuable comments which have indeed made me learnt a lot. I have indeed dedicated three days refining my code however, it appears that tolerance do play a role. I tried the R script which was written as below
if(DEBUG_pH) { #this is on! it calculates pH for each time step
#---------------------------------------
#pH moduleif(ti==2) { y10f <- 0 y20f <- 0 y30f <- 0 y40f <- 0 x10f <- 0 x20f <- 0 x30f <- 0 } InCf <- TANproduction[ti] * (ASL/WSL) /14/2 OutNf <- (SoilEmissionFlux[ti-1] * (PatchArea / 100))/14 + (((TANBudget[ti-1]) * (1-LOSS))/14) InNf <- TANproduction[ti] * (ASL/WSL) /14 Tfix <- T[ti]+273.15 #K Volumefix <- WaterBudget2cm[ti] * (ASL/WSL) #dm3 Volumegfix <- (Porosity * PatchArea * ASL) - Volumefix #dm3 OutCf <- 0 zH0f <- 10^(-pH[ti-1])
Conc_calc ←
function(zHV){
T=Tfix
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=OutCfA <- x10 + x20 + x30 + InN - OutN B <- y10 + y20 + y30 + y40 + InC Kw <- 10^(-14) K1 <- 5.67 * (10 ^ (-10)) * exp(-6286 * ((1 / T) - (1 / 298.15))) K2 <- (10 ^ (-1 * (((2902.39 / T) + (0.02379 * T) - 6.4980)))) K3 <- (10 ^ (-1 * (((3404.71 / T) + (0.032786 * T) - 14.8435)))) H1 <- (56 * exp(4092 * ((1 / T) - (1 / 298.15)))) / ((1 / 8.314) * 0.001 * (1 / T) * (1.013 * 10 ^ 5)) H2 <- ((1/29.41) * exp(2400 * ((1 / T) - (1 / 298.15)))) / ((1 / 8.314) * 0.001 * (1 / T) * (1.013 * 10 ^ 5)) zH <- zHV * V vec1 <- c(0, 0, 0, 0, 0, 1, 0, -zH, 0, 0, (1+H1*(V/Vg)), 0 ) #x1 vec2 <- c(0, 0, 0, 0, 1, 0, 0, K1*V, 0, 0, H1*(V/Vg), 0 ) #x2 vec6 <- c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0) #x3 vec3 <- c(0, 0, 1, 0, 0, 0, 0, 0, -zH, 0, 0, H2*(V/Vg) ) #y1 vec4 <- c(0, 1, 0, 0, 0, 0, 0, 0, K2*V, -zH, 0, H2*(V/Vg) ) #y2 vec5 <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, K3*V, 0, (1+H2*(V/Vg)) ) #y3 vec7 <- c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0) #y4 vec8 <- c(1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0) #RH2 vec9 <- c(-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #RK3 vec10 <- c(0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0) #RK2 vec11 <- c(0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0) #RK1 vec12 <- c(0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0) #RH1 Amatrix <- cbind(vec1, vec2, vec6, vec3, vec4, vec5, vec7, vec8, vec9, vec10, vec11, vec12) colnames(Amatrix)<- c("x1", "x2", "x3", "y1", "y2", "y3", "y4", "RH2", "RK3", "RK2", "RK1", "RH1") Bvector <- c(y30, y20+InC, y10, y40-OutC, x20+InN, x10, x30-OutN, 0, 0, 0, H1*A*(V/Vg), H2*B*(V/Vg)) #const solve(Amatrix,Bvector) roots <- solve(Amatrix,Bvector) zH0 <- zHV0*V pH0 <- -log10(zH0/V) pH <- -log10(zH/V) (zH0 -InC + (-roots[9]+roots[10]+roots[11]) + (buff * (pH-pH0)))-((zH)) }
Blockquote
I see that when I print zHV in the R code and I didn’t specify the tolerance I see that the value of zHV changes. The author of the R code added this in the R code (related to tolerance i.e., tol in R)
Hconc[ti] ← uniroot(Conc_calc, interval=c(0,1), tol=10^-25)$root
When I set the tolerance in R I see that the result looks as below
[1] 0
[1] 1
[1] 0.5
[1] 0.25
[1] 0.125
[1] 0.0625
[1] 0.03125
[1] 0.015625
[1] 0.0078125
[1] 0.00390625
[1] 0.001953125
[1] 0.0009765625
[1] 0.0004882812
[1] 0.0002441406
[1] 0.0001220703
[1] 6.103516e-05
[1] 3.051758e-05
[1] 1.525879e-05
[1] 7.629395e-06
[1] 9.983235e-06
[1] 9.500752e-06
[1] 9.447661e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
> [1] pH = 5.024671 (i.e., -log10 H+ concentration) !which is the default result from the R code
However, in the same R code if the tolerance is not set in the R code the result is as follows
Hconc[ti] ← -uniroot(Conc_calc, interval=c(0,1))$root
pH[ti] ← -log10(Hconc[ti])
if(ti == 2 & ti3 == 1){
print(uniroot(Conc_calc, interval=c(0,1))$root)
print(pH)
[1] 0
[1] 1
[1] 0.5
[1] 0.25
[1] 0.125
[1] 0.0625
[1] 0.03125
[1] 0.015625
[1] 0.0078125
[1] 0.00390625
[1] 0.001953125
[1] 0.0009765625
[1] 0.0004882812
[1] 0.0002441406
[1] 0.0001220703
[1] 0.0001220703
[1] 0.0001220703
[1] 3.91339 (i.e., -log10 H+ concentration)
I have corrected the code implementing major comments for now and printed out variables out here is the updated code in Fortran which there was no setting of tol like in R: and what i printed.
subroutine pH_calculation
use root_module, only: root_scalar
use iso_fortran_env, only: wp => real64
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(wp), 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.0_wp y20f = 0.0_wp y30f = 0.0_wp y40f = 0.0_wp x10f = 0.0_wp x20f = 0.0_wp x30f = 0.0_wp end if InCf = TANproduction(ti, ti3) * (ASL/WSL) / 14.0_wp / 2.0_wp OutNf = (SoilEmissionFlux(ti-1, ti3) * (PatchArea_L/100.0_wp)) / 14.0_wp + (((TANBudget(ti-1, ti3)) * & (1.0_wp-LOSS)) / 14.0_wp) InNf = TANproduction(ti, ti3) * (ASL/WSL) / 14.0_wp Tfix = Ts(ti) + 273.15_wp ! K Volumefix = WaterBudget2cm(ti, ti3) * (ASL/WSL) ! dm3 Volumegfix = (Porosity * PatchArea_L * ASL) - Volumefix ! dm3 OutCf = 0.0_wp zH0f = 10.0_wp ** (-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('bisection',Conc_calc, 0.0_wp, 1.0_wp, zHV, result, iflag) 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) ! the root of the function is the H+ concentration which the author of the R code used !MOA called it roots 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_wp ** (-pH(ti,ti3)) end if end do end do contains function Conc_calc(zHV) result(result) use iso_fortran_env, only: wp => real64 IMPLICIT NONE 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 real(wp), parameter :: TREF=298.15_wp integer :: n=12 ! specify the size of the matrix A = x10 + x20 + x30 + InN - OutN B = y10 + y20 + y30 + y40 + InC Kw = 10.0_wp**(-14.0_wp) ! Kw = 10.0**(-14.0) at 25 C !Print *, 'zHV', Tf, ti K1 = 5.67_wp * (10.0_wp**(-10.0_wp)) * exp(-6286.0_wp * ((1.0_wp / Tf) - (1.0_wp / TREF))) !defines the first H2CO3 dissociation into HCO3- K2 = (10.0_wp**(-1.0_wp * (((2902.39_wp/ Tf) + (0.02379_wp * Tf) - 6.4980_wp)))) !defines the second HCO3- dissociation into CO3-- K3 = (10.0_wp**(-1.0_wp * (((3404.71_wp / Tf) + (0.032786_wp * Tf) - 14.8435_wp)))) !defines the third H2CO3-- dissociation into H2O and co2 H1 = (56.0_wp * exp(4092.0_wp * ((1.0_wp/ Tf) - (1.0_wp / TREF)))) / ((1.0_wp / 8.314_wp) & * 0.001_wp * (1.0_wp / Tf) * (1.013_wp * 10.0_wp ** 5.0_wp)) H2 = ((1.0_wp / 29.41_wp) * exp(2400.0_wp * ((1.0_wp / Tf) - (1.0_wp / TREF)))) & / ((1.0_wp / 8.314_wp) * 0.001_wp * (1.0_wp / Tf) * (1.013_wp * 10.0_wp ** 5.0_wp)) 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 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!There is another function here
end subroutine
If I print the root from the fortran code here is what I obtain. could you please give a feedback.
PS C:\Users\GAG> .\my_program.exe
root is 0.0000000000000000
root is 1.0000000000000000
root is 0.50000000000000000
root is 0.25000000000000000
root is 0.12500000000000000
root is 6.2500000000000000E-002
root is 3.1250000000000000E-002
root is 1.5625000000000000E-002
root is 7.8125000000000000E-003
root is 3.9062500000000000E-003
root is 1.9531250000000000E-003
root is 9.7656250000000000E-004
root is 4.8828125000000000E-004
root is 2.4414062500000000E-004
root is 1.2207031250000000E-004
root is 6.1035156250000000E-005
root is 3.0517578125000000E-005
root is 1.5258789062500000E-005
root is 7.6293945312500000E-006
root is 3.8146972656250000E-006
root is 1.9073486328125000E-006
root is 9.5367431640625000E-007
root is 4.7683715820312500E-007
root is 2.3841857910156250E-007
root is 1.1920928955078125E-007
root is 5.9604644775390625E-008
root is 2.9802322387695312E-008
root is 1.4901161193847656E-008
root is 7.4505805969238281E-009
root is 3.7252902984619141E-009
root is 1.8626451492309570E-009
root is 9.3132257461547852E-010
root is 4.6566128730773926E-010
root is 2.3283064365386963E-010
root is 3.4924596548080444E-010
root is 2.9103830456733704E-010
root is 3.2014213502407074E-010
root is 3.3469405025243759E-010
root is 3.2741809263825417E-010
root is 3.2378011383116245E-010
root is 3.2196112442761660E-010
root is 3.2287061912938952E-010
pH is 9.4909714735041533
I think form the R code it is obvious the tolerance matter especially for a sensitive variable I am simulating? Please could you all give a feedback from my code snippet for the subroutine especially with the root_fortran module? and also any possible solution.
Thanks a lot for your help.
Blockquote