Code Comments

Programming Forum and web based access to our favorite programming groups.
For Programmers: Free Programming Magazines | New: Database administration forum
Registration is free! Edit your profileCalendarFind other membersFrequently Asked QuestionsSearch -> 
Post New Thread











Thread
Author

what is wrong with this code?can someone walk me through this?
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,St
okes
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
Nu = 0.0000181
MP = RHO * ((4.0/3.0)*PI*RP**3.0)
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)
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



Report this thread to moderator Post Follow-up to this message
Old Post
madhu0319@gmail.com
03-27-08 03:53 AM


Re: what is wrong with this code?can someone walk me through this?
In article <2536c6bf-4017-43c3-aa00-5d9eb8dd305e@u72g2000hsf.googlegroups.co
m>,
madhu0319@gmail.com writes:
> Heres a fortran code to solve
....
> 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,St
okes
>                                 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
>                                 Nu = 0.0000181

The problems in this code are too numerous to expend
the effort until you actual post code that has a just
to compile with any compiler.

--
steve

Report this thread to moderator Post Follow-up to this message
Old Post
Steven G. Kargl
03-27-08 03:53 AM


Re: what is wrong with this code?can someone walk me through this?
madhu0319@gmail.com 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,St
okes
>                                 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.

Dick 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
>
>

Report this thread to moderator Post Follow-up to this message
Old Post
Dick Hendrickson
03-27-08 03:53 AM


Re: what is wrong with this code?can someone walk me through this?
On Mar 26, 8:57=A0pm, madhu0...@gmail.com wrote:
> Heres a fortran code to solve the equation of =A0motion 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
> =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 REAL*8 DER=[/color
]
IV_FX,
> DERIV_FY,N,RP,RF,RHO,UF,alpha,K,R0,U
> =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 REAL*8
>  Nu,x,y,y0,THETA,DT,T0,UT0,PI,T,THETA0,St
okes
> =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 REAL*8
>  UR0,b,x1,VX,VY,R,UT_K,UR_K,THETA2,MP,R2,
x_i
> =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 INTEGER t_=[/color
]
max, time
[snip rest of listing]

1. Is this fixed format or free format source code? Your compiler may
truncate lines beyond 72 or 132 characters without warning.
2. Have you used an external diagnostic program to examine your code
like FTNCHEK?
3. Have you used your compiler's diagnostic options?

- e


Report this thread to moderator Post Follow-up to this message
Old Post
e p chandler
03-27-08 03:53 AM


Re: what is wrong with this code?can someone walk me through this?
On Mar 26, 8:19 pm, ka...@troutmask.apl.washington.edu (Steven G.
Kargl) wrote:
> In article <2536c6bf-4017-43c3-aa00-5d9eb8dd3...@u72g2000hsf.googlegroups.
com>,
>         madhu0...@gmail.com writes:
>
>
> 
> .... 
> 
>
> The problems in this code are too numerous to expend
> the effort until you actual post code that has a just
> to compile with any compiler.
>
> --
> steve

Sorry, I did not really understand what u said. I do have an
executable file for this code that works, but i dont know how to post
it here..What do I do abt this program, is it fixable at all?

Report this thread to moderator Post Follow-up to this message
Old Post
madhu0319@gmail.com
03-27-08 03:53 AM


Re: what is wrong with this code?can someone walk me through this?
Thankyou, is there any way i can attach files here? I have an exe file
for the code above, but i dont know if it returns correct values, i
dont anything abt coding. i am a zero. this is another guys work that
i have to improve. i have to get better results and submit it..is this
program feasible in a time frame of 20days?or is it entirely wrong.

heres what is going on in the program..

What basically we are doing is tracking a particles trajectory and
checking if it=A2s hitting a fiber filter that we have. It=A2s a low
Reynolds number flow. That program is not working at very small stokes
number and at every particle velocity. So I have to alter the
program.

For the particle trajectory we are solving an equation

dv/dt=3D 18 * =E7*(U-V) / (2rp )2 =F1p
----------------------------------------    equation 1
where =E7 =3D (2*RP)^2*RHO*UF/ (6*PI*Nu*2*RF)
This is also the equation 4.19 of the textbook that the program
refers to while solving.

The program goes like this

1.	Reading Initialization parameters namely maximum number of
iterations, packing density, particle density, particle radius, fiber
radius and free stream fluid velocity.
2.	Estimate the release point using following relationship: x=3D (b2-
y02) ^0.5 and start iterations.
Here y0 =3D ni*b where b is the radius of the imaginary fluid
surrounding the fiber and ni is a random number designating the
initial position of the particle and 0< ni  < 1.
B=3D RF*alpha^ (-1/2)

3.	Start with ni=3D 1 and keep decreasing it by 1% till the following
condition is satisfied
4.	condition x^2 + y^2 <=3D (rf+rp)^2. Subscript f refers to fiber and
Subscript p refers to particle.  Here x and y are the coordinates
obtained during the execution of the Runge Kutta 4 method for equation
1.

5.	Release particle from the kuwabara cell.
6.	Print output when the particle touches the fiber.
For each iteration we solve equation 1 for X and Y coordinates points
and the velocity components using the RK 4 method. Print x and y
coordinates the velocity components after each iteration.

7.	The coordinates are printed when the particle touches the fiber.

8.	RK 4 method: This is used to calculate the value of of say a
variable say y. By using its present value say yn, we can calculate
the next value of y that is y (n+1).



In our present program since we have to calculate x, y, VX, VY, T:
We have 4 sets of k values by which I mean Our k1, k2, k3, k4 will be.
So we have to solve that equation for x, y , VX. VY, . T is calculated
by simple increment of time step that is T=3DT+DT.


D1x	D1y	D1u	D1v
D2x	D2y	D2u	D2v
D3x	D3y	D3u	D3v
D4x	D4y	D4u	D4v

D1X actually means change in position in a certain time. Its  Dx=3D
Dt*VX ( change in time multiplied by the velocity X which means change
in position ) . U can look at it this way : DX/DT =3D VX (distance in X
direction over time is X velocity). Same is the case with D1Y. D1Y=3D
Dt*VY.

D1U is the change in velocity in the X direction which will be
acceleration. We know that F=3DM*a. Which means acceleration a=3D change
in velocity =3D F/M. Implies change in velocity in the Xdirection is
(forces in Xdirection over M). So D1U/DT=3D FX or
D1U=3D DT*FX. Since F is a sum of all forces F is defined as a function
of x,y, u, v etc

Calculate FX, Fy in separate functions to make things easier . Another
reason  FX and Fy are calculated separate is that they represent the
RHS of equation of motion.




Report this thread to moderator Post Follow-up to this message
Old Post
madhu0319@gmail.com
03-27-08 03:53 AM


Re: what is wrong with this code?can someone walk me through this?
i did not use any diagnostic program. i dont even know anything abt
it. i dont know anything abt fortran or coding..this is my project at
school thats due in 20days..cani run this code on any compiler, i mean
is it compatible with any? shud it return the same values on any
compiler?

pls bear with my silly questions, i am very stressed out abt this..

Report this thread to moderator Post Follow-up to this message
Old Post
madhu0319@gmail.com
03-27-08 03:53 AM


Re: what is wrong with this code?can someone walk me through this?
In article <4e374a6d-1203-40a2-8b15-216aa141ea8d@x41g2000hsb.googlegroups.co
m>,
madhu0319@gmail.com writes:
> On Mar 26, 8:19 pm, ka...@troutmask.apl.washington.edu (Steven G.
> Kargl) wrote: 
>
> Sorry, I did not really understand what u said. I do have an
> executable file for this code that works, but i dont know how to post
> it here..What do I do abt this program, is it fixable at all?

I'm not interested in your executable.

The code you posted isn't in a proper form, so no Fortran compiler
can compile it.  You need to post Fortran code that is in a *readable*
and standard conforming form.

I suggest that you get a few books from your local library on
Fortran programming.  You have 20 days to read the material,
and perhaps learn something.

--
Steve

Report this thread to moderator Post Follow-up to this message
Old Post
Steven G. Kargl
03-27-08 03:53 AM


Re: what is wrong with this code?can someone walk me through this?
On Mar 26, 10:08=A0pm, madhu0...@gmail.com wrote:
> i did not use any diagnostic program. i dont even know anything abt
> it. i dont know anything abt fortran or coding..this is my project at
> school thats due in 20days..cani run this code on any compiler, i mean
> is it compatible with any? shud it return the same values on any
> compiler?
>
> pls bear with my silly questions, i am very stressed out abt this..

You CAN get your program to compile if you

1. remove leading spaces on each line
2. change

EXTERNAL FX, FY, DERIV_FX, DERIV_FY

to

EXTERNAL FX, FY

3. compile as free format source

DERIV_FX  and DERIV_FY are variables NOT functions passed as arguments
to routines so there is no need to make them EXTERNAL.

Making the numerical changes suggested by others DOES not make any
difference in the output except for trailing digits for the value of
Stokes.

When I compiled and ran the program it printed:

*** Parameters read from file particle.par
***
The Stokes is equal to:
0.0006589939662307433
n_t=3D 1. y=3D 0.000023255106965997813
n_t=3D 0.995 y=3D 0.000023138831431167825
n_t=3D 0.990025 y=3D 0.000023023137274011985
n_t=3D 0.985074875 y=3D 0.000022908021587641927
n_t=3D 0.9801495006250001 y=3D 0.000022793481479703716
n_t=3D 0.9752487531218751 y=3D 0.0000226795140723052
n_t=3D 0.9703725093562657 y=3D 0.000022566116501943672
n_t=3D 0.9655206468094844 y=3D 0.000022453285919433955
n_t=3D 0.960693043575437 y=3D 0.000022341019489836785

Note that your value of Stokes is NOT between 0.1 and 2.0 as desired.
Maybe better input data will help. Maybe better values for other
constants will help.

Check the program listing against your formulas and equations. Verify
that your program is correctly transcribed.
Insert WRITE statements in the program to display the value of
variables as the program runs. Compare these to manual calculations.

What are reasonable ranges for the input variables? Try combinations
of LOW, MID and HIGH values for the input variables and see what you
get.

At this point there is no substitute for pouring over your program
listing and output with paper, pencil and a calculator.

With MINOR changes the program DOES compile. You don't have obvious
errors involving mismatched arguments to subroutines or common blocks
that do not match. You may have other

programming errors
logical errors
bad numerical behavior.

Other than the "Non-Fortran" flavor of your program, there is nothing
obviously wrong that jumps out except

goto 5
followed by
STOP. :-).

Good luck.

-- e


Report this thread to moderator Post Follow-up to this message
Old Post
e p chandler
03-27-08 10:14 AM


Re: what is wrong with this code?can someone walk me through this?
If the program becomes exhausted and can't make it to 5, maybe it will
just fall into STOP.

Good suggestions.  Reminds me of my old old old old college days.

On Wed, 26 Mar 2008 21:21:53 -0700 (PDT), e p chandler <epc8@juno.com>
wrote:
>Other than the "Non-Fortran" flavor of your program, there is nothing
>obviously wrong that jumps out except
>
>goto 5
>followed by
>STOP. :-).

Report this thread to moderator Post Follow-up to this message
Old Post
Richard L Walker
03-27-08 10:14 AM


Sponsored Links




Last Thread Next Thread Next
Pages (3): [1] 2 3 »
Search this forum -> 
Post New Thread

Fortran archive

Show a Printable Version Send to friend Email This Page to Someone! subscribe to this thread Receive updates to this thread
Computer Consultants
Programming Jobs
Visual Basic Controls
SQL Server Programming
Webservices
Java Security
Visual Studio
C# Programming
Visual J++
Software engineering
Open source Software
Perl Programming
PHP Programming
ASP Programming
ASP .NET Programming
Visual Basic Programming
Windows Scripting Host
Java Programming
Java Help
Java Beans
VBScript
Cobol
MAC Applications
Unix Programming
Forum Jump:
All times are GMT. The time now is 02:11 PM.

 
Free MCSE Braindumps | Real Estate Topics

Programming forum archive

Copyrights CodeComments.com 2004 - 2006

Powered by vBulletin Copyright 2000-2006 Jelsoft Enterprises Limited.