0

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
Jake
  • 13
  • 5

2 Answers2

1

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.

chw21
  • 7,970
  • 1
  • 16
  • 31
  • Ah, so I need to insert those other arguments into each function call statement? – Jake Dec 14 '15 at 00:22
  • You probably mean either to call the function `fx` or the result `f`, depending on which interface you are copying. It's possibly also worth noting the mismatch between the `intent(inout)` and `intent(out)` arguments in terms of which is defined when it comes to the call to `secant`. [And that the real named constants don't really belong in the interface block.] – francescalus Dec 14 '15 at 00:27
  • Yes I changed the secant function calls to something like fx(xold,L,Q,hf,rho,mu,rough) and im getting no implicit type errors. I'm not exactly sure what you mean regarding the INOUT. The xold and xolder are being passed in from the main program, and the xnew will be passed back as the answer. Should I change the INOUT to just IN? – Jake Dec 14 '15 at 00:37
  • @Jake Sorry for confusion. The first part of my comment was pointing out a potential typo in this answer. The second part, about the intents, was a suggestion for improvement: `xolder` is associated with an `intent(inout)` argument without first being defined, and `xnew` is defined but associated with an `intent(out)` argument. There is perhaps a mismatch there. – francescalus Dec 14 '15 at 00:52
  • @francescalus thank you, after trying to fix that I just moved all declaration statements and variable read ins to the secant module and it compiled. Now I am getting negative answers for diameter so I have to take a look at where the error is occurring. I'll update the code to see if anyone could point out any new mistakes. – Jake Dec 14 '15 at 01:34
  • @Jake, it's probably best not to change that question looking for more errors after using the suggestion here. If this answer solves the problem, at least as you stated it first, you should accept it. You can ask another question for subsequent points. – francescalus Dec 14 '15 at 01:39
0

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.

roygvib
  • 7,218
  • 2
  • 19
  • 36
  • thank you for taking the time out of your day to help me. I originally solved this problem in excel using solver and was given 0.799 inches, reasonable for the size of this system (according to my mentor). As far as the solution goes, I've tried cleaning up the parenthesis and ended up getting the same answers. Units have also been checked multiple times, and I am certain that those are right. I rewrote the function using my empirical derivation which took all factors besides D out of the function and got the same thing. I didn't want to hard-wire the program and limit its usefulness. – Jake Dec 14 '15 at 04:44
  • @Jake Hmm, I see... and if you can calculate the answer by Excel, does it agree with the Fortran solution above (e.g. after some unit conversion etc)? At least, the secant code itself seems to be working correctly. [Also, by "pre-calculate the constant factors", I mean using some temporary variables to store constant factors as a whole, e.g. `coeffA = hf/(L* 4*Q/pi)` to make it easy to see the structure of the expression (not hard-coding the literal numbers). This is also useful for the reader of this site, because it takes time to "decode" the program... – roygvib Dec 14 '15 at 05:25
  • 1
    Yes on excel the answer is reported in feet, so the number is multiplied by 12 to get to inches which ends up being ~0.799. I thought about doing coefficients but thought I couldn't because it would mess the function up but I am working on that now. I also realized I forgot a few terms from the colebrook equation which would explain a lot I am currently working with f=(1/SQRT(hf/((L/D)*(((4.0*Q)/(pi*D**2))/2*g)))) & +2.0*log10((rough/(3.7*D))+ & ((2.51/((rho*(4.0*Q/pi*D**2)*D/mu) & *SQRT(hf/((L/D)*(((4.0*Q)/(pi*D**2))/2*g))))))) But am going to simplify it. Thanks for your help – Jake Dec 14 '15 at 05:30
  • Well some exponent stars didnt go through, but I was forgetting the 2g term on my velocity head and the fact that flow rate / area gives a D^2 value in the denominator. – Jake Dec 14 '15 at 05:33
  • @Jake Oh, if you want to write equations, it would be definitely better to modify the original question (not in the comments), because the typeset is very poor here... But you can enclose an equation with backquote (`) to highlight it. Also plz note that Fortran has log() and log10() (the base is different), and also it is more accurate to use double-precision literal constants like 2.51d0. – roygvib Dec 14 '15 at 05:36
  • Also I was confused as to the listing of formal arguments if I use coefficients. My coefficients I'm going to create will be friction factor, and reynolds number, which use most of the listed variables. Will I have to change the argument list of the function to reflect this, since these variables are being used directly in the function calculation? – Jake Dec 14 '15 at 05:39
  • Usually, it is simplest to pass the parameters (like Q, L, etc) to the function by defining a `module`, placing those parameters in that module, and `use`-ing them from any subroutine that needs them. I remember there are several good Q&A in this site about how to do this (please use the search window right above, e.g., "fortran module". – roygvib Dec 14 '15 at 05:49
  • I've updated the code, trying smaller guesses now. I originally have shared variables in a module but was getting errors related to function definitions and decided to just keep it within modules. I've rewritten this thing more than a few times haha. Thank you for all your help, I may just give up on this and turn in my report as is. I was supposed to study for a final today and got stuck fixing this for the last 10 hours. – Jake Dec 14 '15 at 06:08
  • Hmm, 10 hours are tough... With your updated equation, the program converges to 0.3794339, so much closer to Excel's result (I used xold=1.0 and xolder=0.5 as starting values). – roygvib Dec 14 '15 at 06:18
  • Hey thanks for verifying that, it will give me confidence in my reasoning for not using this answer. Doesn't seem to change within reasonable guess limits. – Jake Dec 14 '15 at 06:40