Re: what is wrong with this code?can someone walk me through this?

madhu0319@xxxxxxxxx wrote:
Heres a fortran code to solve the equation of motion of very small
particle in a low reynolds number flow. It is not working in all
compilers, I have a working version of this code, but it doesnt work
for all inputs. Can someone tell me what is wrong with this and fix it
so that it works for a stokes number of 0.1 through 2. Stokes number
is calculated in the program. You need to work backwards by playing
with the input file and see for what values of input the stokes number
is 0.1 through 2.

my particle hardly moves. I have no way of even telling if my values
are correct. What should I do? Pls help.
my input file. it contains values of t_max, alpha, Rho, Rp, Rf, Uf in
that order exactly
example working input file values are:
600,0.05,2165,3.00E-07,5.20E-06,0.003 in that order.

source code:

IMPLICIT NONE
REAL*8 DERIV_FX,
DERIV_FY,N,RP,RF,RHO,UF,alpha,K,R0,U
REAL*8
Nu,x,y,y0,THETA,DT,T0,UT0,PI,T,THETA0,Stokes
REAL*8
UR0,b,x1,VX,VY,R,UT_K,UR_K,THETA2,MP,R2,x_i
INTEGER t_max, time
REAL*8 n_t
EXTERNAL FX, FY, DERIV_FX, DERIV_FY
COMMON Nu, RP, RF, MP, UF, alpha, K,
PI, Stokes, x_i
CALL read_parameters (t_max, alpha,
RHO, RP, RF, UF)
PI = 3.1415
I don't know enough about your problem to debug it. But, here's
a simple problem. You've declared all of you variables to be
real*8, which is usually called double precision. But, you specify
all of your constants as single precision. 3.14159 is a single
precision value. The fact that it is assigned to a double precision
variable does not make it more precise; values are usually converted
from single to double by appending zeros in the low order bits. It
doesn't make much difference for PI, since you only specify a few
digits.

For any number that isn't an integer or a "power of 2" fraction
like .5 or .25, specify the value in double precision. The way
to do that, consistent with your coding style, is to write
3.1415D0
Better yet to look up a few more digits for PI, it can't hurt
to be more accurate ;) .

Nu = 0.0000181
MP = RHO * ((4.0/3.0)*PI*RP**3.0)
This is a potentially worse problem. The expression (4.0/3.0)
will be computed as approximately 1.333333 and then converted to
double precision and become approximately 1.3333330000000. You
probably wanted 1.3333333333333.

For expressions like (4.0/3.0), make sure at least one of the values
is double precision. Best to do (4.0D0/3.0D0) .

Secondly, you almost never want to use expressions like
RP**3.0 where the exponent is a real number. Most compilers
will do that with log and exp functions and lose accuracy.
Better to do exponentiation as RP**3 when the exponent is an
integer.

Stokes = (2.0*RP) **2*RHO*UF/
(6*PI*Nu*2.0*RF)
write (6,*) 'The Stokes is equal to:'
write (6,*) Stokes
b = RF * alpha**(-1.0/2.0)
Personally, I'd use the SQRT function, rather than exponentiation to .5.
It's likely to be faster, more accurate, and much easier for humans
to read.

*** Hendrickson
K =
-0.5*log(alpha)-0.75+alpha-0.25*alpha**2.0
n_t = 1.0
5 y0 = n_t*b
x= (-1.0)*(b**2.0 -
y0**2.0)**(1.0/2.0)
R0 = (x**2.0 + y0**2.0)**(1.0/2.0)
THETA0 = atan(y0/x)


THETA0= PI/2.0 + THETA0
DT = 0.000000001
T0 = 0.0
T =T0
THETA = THETA0
R = R0
Y= Y0
write (6,*) 'n_t=', n_t, 'y=', y
U = UF

UT0 = (U/(2.0*K))*(2.0*log(R/RF)+1+alpha-(RF**2.0/R**2.0) *(1.0-alpha/
2.0)-(3.0*alpha/2.0)*(R**2.0/RF**2.0))*sin(THETA)
UR0 = (U/(2.0*K))*(2.0*log(R/RF)-1.0+alpha+(RF**2.0/R**2.0) *(1.0-
alpha/2.0)-(alpha/2.0)*(R**2.0/RF**2.0))*cos(THETA)
VX= UT0
VY = UR0
OPEN (11, FILE='PARTICLEOUT.DAT')
WRITE (11,*) 'OUTPUT VARIABLES =
time,THETA, R, x, y, VT, UT,VR, UR'
DO 7 time = 1, t_max
x_i =x
CALL KUTTA(x,y,VX,VY,T,DT,FX,FY)
THETA2 = atan(y/x)
R2 = (x**2.0 + y**2.0)**(1.0/2.0)
UT_K = (U/(2.0*K))*(2.0*log(R2/RF)+1.0+alpha-(RF**2.0/R2**2.0)*(1.0-
alpha/2.0)-3.0*alpha/2.0*(R2**2.0/RF**2.0))*sin(THETA2)
UR_K = (U/(2.0*K))*(2.0*log(R2/RF)-1.0+alpha+(RF**2.0/R2**2.0)*(1.0-
alpha/2.0)-alpha/2.0*(R2**2.0/RF**2.0))*cos(THETA2)
WRITE(11,10)
time,THETA2,R2,x,y,VX,UT_K,VY,UR_K

10 FORMAT(I1,9E15.6)

if ((x**2+y**2) .LE. (RF+RP)**2) then
write (6,*) 'The Stokes is equal to:'
write (6,*) Stokes
stop
endif

if (n_t .LE. 8.0E-8) then
close (11)
write (6,*) 'Neglegible single fiber
efficiency under these conditions'

stop
endif

if ((x-x_i) .LE. 0.0) then
n_t = n_t - (n_t/200.0)
close (11)
goto 5
stop
endif
7 CONTINUE
CLOSE (11)
END

SUBROUTINE KUTTA(x,y,VX,VY,T,DT,FX,FY)
IMPLICIT NONE
REAL*8
D1X,D1Y,DT,VX,VY,D1U,D1V,D2X,D2Y,D2U,D2V
REAL*8
D3X,D3Y,D3U,D3V,D4X,D4Y,D4U,D4V,T,x,y
REAL*8 DERIV_FX, DERIV_FY

DERIV_FX = 0.0
DERIV_FY =0.0

D1X = DT * VX
D1Y = DT * VY

CALL FX(x,y,VX,DERIV_FX)
CALL FY(x,y,VY,DERIV_FY)

D1U = DT * DERIV_FX
D1V= DT * DERIV_FY

D2X = DT * (VX+D1U/2.0)
D2Y = DT * (VY+D1V/2.0)

CALL FX(x+D1X/2.0, y+D1Y/2.0, VX+D1U/
2.0, DERIV_FX)
CALL FY(x+D1X/2.0, y+D1Y/2.0, VY+D1V/
2.0, DERIV_FY)

D2U = DT * DERIV_FX
D2V = DT * DERIV_FY

D3X = DT * (VX+D2U/2.0)
D3Y = DT * (VY+D2V/2.0)

CALL FX(x+D2X/2.0, y+D2Y/2.0, VX+D2U/
2.0, DERIV_FX)
CALL FY(x+D2X/2.0, y+D2Y/2.0, VY+D2V/
2.0, DERIV_FY)

D3U = DT * DERIV_FX
D3V = DT * DERIV_FY

D4X = DT * (VX+D3U)
D4Y = DT * (VY+D3V)

CALL FX(x+D3X, y+D3Y, Vx+D3U,
DERIV_FX)
CALL FY(x+D3X, y+D3Y, VY+D3V,
DERIV_FY)

D4U = DT * DERIV_FX
D4V= DT * DERIV_FY

T = T + DT
X= x + (D1X + 2.0*D2X + 2.0*D3X +
D4X) / 6.0
y =y + (D1Y + 2.0*D2Y + 2.0*D3Y +
D4Y) / 6.0
VX = VX + (D1U + 2.0*D2U + 2.0*D3U +
D4U) / 6.0
VY = VY + (D1V + 2.0*D2V + 2.0*D3V +
D4V) / 6.0

RETURN
END

SUBROUTINE FX (x, y, VX, DERIV_FX)

IMPLICIT NONE
REAL*8
UT_K,Sk,K,Nu,alpha,y,RP,RF,MP,UF,x,VX
REAL*8
Stokes,DERIV_FX,PI,THETA,R,x_i,U
COMMON
Nu,RP,RF,MP,UF,alpha,K,PI,Stokes,x_i

R= (x**2.0 + y**2.0)**(1.0/2.0)
U = UF

THETA = atan(y/x)
UT_K = (U/(2.0*K))*(2.0*log(R/RF)
+1.0+alpha-(RF**2.0/R**2.0)*(1.0-alpha/2.0)-(3.0*alpha/2.0)*(R**2.0/
RF**2.0))*sin(THETA)
Sk = 6.0*PI*Nu*RP/MP

DERIV_FX = Sk*(UT_K-VX)

RETURN
END

SUBROUTINE FY(x, y, VY, DERIV_FY)

IMPLICIT NONE
REAL*8 UR_K, Sk, K, Nu, alpha, y, RP, RF, MP,
UF, THETA, VY
REAL*8 Stokes,DERIV_FY,PI,x,R,x_i,U
COMMON Nu,RP,RF,MP,UF,alpha,K,PI,Stokes,x_i

R= (x**2.0 + y**2.0)**(1.0/2.0)
U = UF


THETA = atan(y/x)
UR_K=( (U/(2.0*K))*(2.0*log(R/RF)-1.0+alpha+
(RF**2.0/R**2.0)*(1.0-alpha/2.0)-(alpha/2.0)*(R**2.0/
RF**2.0))*cos(THETA))

Sk = 6.0*PI*Nu*RP/MP


DERIV_FY = Sk* (UR_K-VY)

RETURN
END

SUBROUTINE read_parameters (t_max,alpha, RHO,
RP, RF, UF)
implicit none
integer t_max
real*8 alpha,RHO,RP,RF,UF
open (10, file='particle.par')

read(10,*) t_max

read(10,*) alpha

read(10,*) RHO

read(10,*) RP

read(10,*) RF

read(10,*) UF

close (10)

write (6,*) '*** Parameters read from file
particle.par'
write (6,*) '***'

return
end


.