Secant method solving for pipe diameter

Question

I am trying to write a program to solve for pipe diameter for a pump system I've designed. I've done this on paper and understand the mechanics of the equations. I would appreciate any guidance.

EDIT: I have updated the code with some suggestions from users, still seeing quick divergence. The guesses in there are way too high. If I figure this out I will update it to working.

MODULE Sec
CONTAINS

SUBROUTINE Secant(fx,xold,xnew,xolder)
IMPLICIT NONE
INTEGER,PARAMETER::DP=selected_real_kind(15)
REAL(DP), PARAMETER:: gamma=62.4
REAL(DP)::z,phead,hf,L,Q,mu,rho,rough,eff,pump,nu,ppow,fric,pres,xnew,xold,xolder,D
INTEGER::I,maxit

INTERFACE
omitted
END INTERFACE

Q=0.0353196
Pres=-3600.0
z=-10.0
L=50.0
mu=0.0000273 
rho=1.940
nu=0.5
rough=0.000005
ppow=412.50
xold=1.0
xolder=0.90
D=11.0
phead = (pres/gamma)
pump = (nu*ppow)/(gamma*Q)
hf = phead + z + pump

maxit=10
I = 1

DO
xnew=xold-((fx(xold,L,Q,hf,rho,mu,rough)*(xold-xolder))/ &
      (fx(xold,L,Q,hf,rho,mu,rough)-fx(xolder,L,Q,hf,rho,mu,rough)))

xolder = xold
xold = xnew
I=I+1
WRITE(*,*) "Diameter = ", xnew
IF (ABS(fx(xnew,L,Q,hf,rho,mu,rough)) <= 1.0d-10) THEN
EXIT
END IF

IF (I >= maxit) THEN
EXIT
END IF 
END DO

RETURN

END SUBROUTINE Secant
END MODULE Sec

PROGRAM Pipes
USE Sec
IMPLICIT NONE
INTEGER,PARAMETER::DP=selected_real_kind(15)
REAL(DP)::xold,xolder,xnew

INTERFACE
omitted
END INTERFACE

CALL Secant(f,xold,xnew,xolder)

END PROGRAM Pipes

FUNCTION f(D,L,Q,hf,rho,mu,rough)
IMPLICIT NONE
INTEGER,PARAMETER::DP=selected_real_kind(15)
REAL(DP), PARAMETER::pi=3.14159265d0, g=9.81d0
REAL(DP), INTENT(IN)::L,Q,rough,rho,mu,hf,D
REAL(DP)::f, fric, reynold, coef

fric=(hf/((L/D)*(((4.0*Q)/(pi*D**2))/2*g)))

reynold=((rho*(4.0*Q/pi*D**2)*D)/mu)

coef=(rough/(3.7d0*D))

f=(1/SQRT(fric))+2.0d0*log10(coef+(2.51d0/(reynold*SQRT(fric))))

END FUNCTION

Answer 1

It seems that the initial values for xold and xolder are too far from the solution. If we change them as

xold   = 3.0d-5
xolder = 9.0d-5

and changing the threshold for convergence more tightly as

IF (ABS(fx(xnew,L,Q,hf,rho,mu,rough)) <= 1.0d-10) THEN

then we get

...
Diameter =    7.8306011049894322E-005
Diameter =    7.4533171406818087E-005
Diameter =    7.2580746283970710E-005
Diameter =    7.2653611474296094E-005
Diameter =    7.2652684750264582E-005
Diameter =    7.2652684291155581E-005

Here, we note that the function f(x) is defined as

FUNCTION f(D,L,Q,hf,rho,mu,rough)
...
f = (1/(hf/((L/D)*((4*Q)/pi*D))))                                   !! (1)
    + 2.0 * log(  (rough/(3.7*D)) + (2.51/(((rho*((4*Q)/pi*D))/mu)  !! (2)
                    * (hf/((L/D)*((4*Q)/pi*D)))))                   !! (3)
               )
END FUNCTION 

where terms in Lines (1) and (3) are both constant, while terms in Line (2) are some constants over D. So, we see that f(D) = c1 - 2.0 * log( D / c2 ), so we can obtain the solution analytically as D = c2 * exp(c1/2.0) = 7.26526809959e-5, which agrees well with the numerical solution above. To get a rough idea of where the solution is, it is useful to plot f(D) as a function of D, e.g. using Gnuplot.

But I am afraid that the expression for f(D) itself (given in the Fortran code) might include some typo due to many parentheses. To avoid such issues, it is always useful to first arrange the expression for f(D) as simplest as possible before making a program. (One TIP is to extract constant factors outside and pre-calculate them.)

Also, for debugging purposes it is sometimes useful to check the consistency of physical dimensions and physical units of various terms. Indeed, if the magnitude of the obtained solution is too large or too small, there might be some problem of conversion factors for physical units, for example.

Answer 2

You very clearly declare the function in the interface (and the implementation) as

FUNCTION f(L,D,Q,hf,rho,mu,rough)
    IMPLICIT NONE
    INTEGER,PARAMETER::DP=selected_real_kind(15)
    REAL(DP), PARAMETER::pi=3.14159265, g=9.81
    REAL(DP), INTENT(IN)::L,Q,rough,rho,mu,hf,D
    REAL(DP)::fx
END FUNCTION

So you need to pass 7 arguments to it. And none of them are optional.

But when you call it, you call it as

xnew=xold-fx(xold)*((xolder-xold)/(fx(xolder)-fx(xold))

supplying a single argument to it. When you try to compile it with gfortran for example, the compiler will complain for not getting any argument for D (the second dummy argument), because it stops with the first error.