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

On Mar 28, 3:09=A0pm, madhu0...@gmail.com wrote:
> On Mar 29, 12:03 am, madhu0...@gmail.com wrote:
>
>
>
>
> 
> 
enews:fsesmg$8a6$1@gnus01.u.washington.edu...
> 
 
> 
L*8 DERIV_FX, 
L*8 
L*8 
EGER t_max, time 
L*8 n_t 
ERNAL FX, FY, DERIV_FX, DERIV_FY 
MON Nu, RP, RF, MP, UF, alpha, K, 
L read_parameters (t_max, alpha, 
=3D 3.1415 
=3D 0.0000181
> 
> 
> 
think 
> 
or the 
> 
> 
> 
> 
> 
> 
> 
>
> Stokes number is defined in the program. The actual definition for
> stokes is, it is ratio of the stopping distance of a particle to a
> characteristic dimension of the obstacle.
>
> In basic terms, its a way of measuring how good the program is..say
> for example instead of saying that the code doesnt work for a certain
> particle radius, we say it doesnt work for this particular stokes
> number and we do this by putting the radius in the formula for stokes
> number, calculate the stokes number (given the radius) and say that it
> is not working (the program) for this particular stokes number..its
> like a standard..in viscous flows
>
> In program terms, i should be getting output of stokes between 0.1 and
> 2.0 for certain values of input but unfortunately i am not getting
> that, am trying to figure out why. [Baffled]- Hide quoted text -
>
> - Show quoted text -

Looking at your program listing, I think I've found the "smoking gun".
Your very long expressions for UT0, UR0, UT_K, UR_K, UT_K and UR_K
again look similar, but have different structures. You've used
parentheses where you don't need them and the rest of your parentheses
get confusing.

Long ago I was taught this trick for figuring out precedence in a
complicated expression with many levels of parentheses.

For each statement

1. start a counter at 0
2. for each left parenthesis add 1 to the counter
3. for each right parenthesis subtract 1 from the counter

Things at the same level are part of the same expression.

Notice that in your code the levels are not consistent between
statements. Often programs have parallel structure. Even though you
may interchange variables and multiply by different constants, the
general form tends to be similar but not so in your program.

Confusing the reader with un-needed REAL exponents and NOT using
SQRT() where it should apply makes this possible problem LESS obvious.

I hope this helps.

--- e

madhu0319@gmail.com wrote:
>
> In program terms, i should be getting output of stokes between 0.1 and
> 2.0 for certain values of input but unfortunately i am not getting
> that, am trying to figure out why. [Baffled]

My approach for dealing with problems like this is to do a parallel
calculation using pencil and paper.  This allows me to break the problem
down into subsets with known answers.  You can modify the program to
display the results of intermediate calculations, and compare with your
p&p results; you may find that simply performing the reduction
identifies the problem for you, otherwise you will likely be able to
localize it in a small piece of the calculation where it's fairly easy
to see what's wrong.

On Mar 31, 6:58=A0pm, Craig Powers <eni...@hal-pc.org> wrote:
> madhu0...@gmail.com wrote:
> 
>
> My approach for dealing with problems like this is to do a parallel
> calculation using pencil and paper. =A0This allows me to break the problem=[/color
]

> down into subsets with known answers. =A0You can modify the program to
> display the results of intermediate calculations, and compare with your
> p&p results; you may find that simply performing the reduction
> identifies the problem for you, otherwise you will likely be able to
> localize it in a small piece of the calculation where it's fairly easy
> to see what's wrong.

One last try at beating this dead horse. :-(.

[to the OP]

A.

1. Is

Stokes =3D (2.0*RP) **2*RHO*UF/ (6*PI*Nu*2.0*RF)

the correct formula?

Not that some obvious algebraic simplification has NOT been done
particularly in the denominator.

Should the 2.0 be an exponent or remain as a factor?

2. Are the units correct?

3. Second above advice to check this formula by hand.

B. One possible reason that the particle does not go anywhere is
fairly obvious. The time step is small =3D 1E-9. Maybe a larger time
step would help. But more bothersome is that the GOTO 5 statement is
executed at all, and in particular for most values of n_t at time step
2.

1. Instrumenting the program with stategic WRITE statements might
help.

C. So far NO ONE has commented on the OP's use of REAL*8 to specify
DOUBLE PRECISION. I guess that the lack of comment indicates that this
usage is the least of the OP's problems here. :-<<<.

-- e [far past bedtime for me...]

"e p chandler" <epc8@juno.com> wrote in message
news:7b2a9762-1136-432d-99ca-74d9a04cbd3e@z38g2000hsc.googlegroups.com...
On Mar 31, 6:58 pm, Craig Powers <eni...@hal-pc.org> wrote:
> madhu0...@gmail.com wrote:
> 
>
> My approach for dealing with problems like this is to do a parallel
> calculation using pencil and paper. This allows me to break the problem
> down into subsets with known answers. You can modify the program to
> display the results of intermediate calculations, and compare with your
> p&p results; you may find that simply performing the reduction
> identifies the problem for you, otherwise you will likely be able to
> localize it in a small piece of the calculation where it's fairly easy
> to see what's wrong.

One last try at beating this dead horse. :-(.

[to the OP]

A.

1. Is

Stokes = (2.0*RP) **2*RHO*UF/ (6*PI*Nu*2.0*RF)

the correct formula?

Not that some obvious algebraic simplification has NOT been done
particularly in the denominator.

Should the 2.0 be an exponent or remain as a factor?

2. Are the units correct?

3. Second above advice to check this formula by hand.

B. One possible reason that the particle does not go anywhere is
fairly obvious. The time step is small = 1E-9. Maybe a larger time
step would help. But more bothersome is that the GOTO 5 statement is
executed at all, and in particular for most values of n_t at time step
2.

1. Instrumenting the program with stategic WRITE statements might
help.

C. So far NO ONE has commented on the OP's use of REAL*8 to specify
DOUBLE PRECISION. I guess that the lack of comment indicates that this
usage is the least of the OP's problems here. :-<<<.

---> OP is lazy.  Fortran can't help him.

I was hoping that we would get to discuss diff eq's rendered in fortran.

Can you bring up this topic again in a month?
--

"That this social order with its pauperism, famines, prisons, gallows,
armies, and wars is necessary to society; that still greater disaster
would ensue if this organization were destroyed; all this is said only
by those who profit by this organization, while those who suffer from it
- and they are ten times as numerous - think and say quite the contrary."

~~  Leo Tolstoy

e p chandler <epc8@juno.com> wrote:

> C. So far NO ONE has commented on the OP's use of REAL*8 to specify
> DOUBLE PRECISION. I guess that the lack of comment indicates that this
> usage is the least of the OP's problems here. :-<<<.

Yeah. Though it is non-standard, it is pretty far down on the list of
things I'd comment about. After first worrying about getting the right
answers, I'd probably comment on some things that actually are standard
conforming, but shouldn't be done anyway...such as pretty much all of
the exponentiations in the forms used; someone already mentioned the
**(-1.0/2.0), but I don't recall if anyone also mentioned things like
**3.0.

But it seemed to me that the OP mostly needed much more fundamental help
along the lines some people have been suggesting (print out intermediate
results and hand-check them).

--
Richard Maine                    | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle           |  -- Mark Twain

On 2008-04-01, Richard Maine <nospam@see.signature> wrote:

> But it seemed to me that the OP mostly needed much more fundamental help
> along the lines some people have been suggesting (print out intermediate
> results and hand-check them).

Another thing that I'd recommended to the OP in private E-Mail is
to check his step sizes.  He'd been using Runge-Kutta with fixed
step sizes and without error estimates, which is a bad idea.

The least he should do is to vary step sizes and determine a
plateau where the step size doesn't influence the results too much.

On Apr 1, 1:29=A0am, nos...@see.signature (Richard Maine) wrote:
> e p chandler <e...@juno.com> wrote:
> 
>
> Yeah. Though it is non-standard, it is pretty far down on the list of
> things I'd comment about. After first worrying about getting the right
> answers, I'd probably comment on some things that actually are standard
> conforming, but shouldn't be done anyway...such as pretty much all of
> the exponentiations in the forms used; someone already mentioned the
> **(-1.0/2.0), but I don't recall if anyone also mentioned things like
> **3.0.
>
> But it seemed to me that the OP mostly needed much more fundamental help
> along the lines some people have been suggesting (print out intermediate
> results and hand-check them).
>
> --
> Richard Maine =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0| Good judgement come=[/color
]
s from experience;
> email: last name at domain . net | experience comes from bad judgement.
> domain: summertriangle =A0 =A0 =A0 =A0 =A0 | =A0-- Mark Twain

In this piece of code the OP has two problems:

n_t =3D 1.0
5  y0 =3D n_t*b
x=3D (-1.0)*SQRT(b**2 - y0**2)
R0 =3D SQRT(x**2 + y0**2)
THETA0 =3D atan(y0/x)
THETA0=3D PI/2.0 + THETA0
=2E..
THETA =3D THETA0
R =3D R0
Y=3D Y0
write (6,*) 'n_t=3D', n_t, 'y=3D', y
U =3D UF
UT0 =3D (U/(2.0*K))*(2.0*log(R/RF)+1+alpha-(RF**2/R**2) *(1.0-alpha/2.0)-
(3.0*alpha/2.0)*(R**2/RF**2))*sin(THETA)
UR0 =3D (U/(2.0*K))*(2.0*log(R/RF)-1.0+alpha+(RF**2/R**2) *(1.0-alpha/
2.0)-(alpha/2.0)*(R**2/RF**2))*cos(THETA)

1. on the first pass, X =3D -0
2. Why calculate sin(pi/2+arctan(y/x)) and cos(pi/2+arctan(y/x))
instead of X/R and -Y/R? (phase shift means take derivatives).

Simple trigonometry avoids the divide by 0 problem, which is made
worse because atan(-0) + pi/2 =3D -pi/2 (double) + pi/2 (single) here.

Later substitutions for sin(THETA) and cos(THETA) are straightforward.

As another poster suggested, decreasing the step size does seem to
help. With DT=3D1E-5, the program still prints junk, but at least the
particle changes (x,y). :-).

-- e