Some help with a Fortran source code

[bigduke88]

Hello everyone
I would like to ask you for some help with a Fortran source code.
I have to make a program, which solves numerically an ordinary differential equation and after that an integral.
The equation, is going to be solved using the finite difference approximation.
This is equation "developed" in finite differences:
Ui+1 = 0.998*h - 0.7366*(Ui)^2*Cd*h + Ui
This is an equation, which describes the velocity of a moving spherical particle.
"Cd" is a drag coefficient and is function of U:
And this is the equation, which describes the drag coefficient:
Cd = (0.087108/U) + (0.622126/U^0.28)

I'm quite new in any use of an algorythmic "language", but I have to solve the equation and then an intergral in Fortran.
Could you help with some advice, how to do this?
Thank you!!!

 

[gummibaer]

And you expect us to deliver the code ??

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[bigduke88]

No, I already have created a code, for some pde-s and an itegral, but if I include this equation in the "original" source code for the pde-s I think, that I'm going to mess things, because I'm working with a dimensionless time, but for the ODE, have to transfer from dimensionless to a "normal" one... I thought about a subroutine - to solve the ode with a different program and then "call" the result in the original one

 

[gummibaer]

Why don't you post what you have ?

I do not have any idea what a PDE or ODE might be (partial differential equation, ordinary differential equation ??) and how this is introduced into your code and what you want to do with it in terms of coding.

With some code we would have a baseline for discussion.

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[salgerman]

As you just realized, in order to program something, you really, REALLY need to know what you are doing.

Contrary to popular believe, computer programs are not magic...other than artificial intelligence, neural networks and other stuff that I don't know...for the rest of us, computers suffer of artificial stupidity and we need to tell them exactly what to do.

But don't let programming scare you...at the minimum, think of programming as a way of automating your hand calculations...so, how would you go about solving this problem by hand? do it, at least once; or, if it is tooo long, sketch the steps up, take a look at the steps and program them.

Lastly, do not expect any of us be experts on the field that you are and so we may not be able to help you go from "dimensionless" time to a "normal" one...if we are experts on something is fortran and whatever other field we use fortran for.

So...think and ask the question in a way we can understand...try to make it a fortran question

 

[fanta2]

Your equation seems to be the final product (explicit solution) after you have substituted the finite difference form into the PDE. If that is the case you have to just calculate Ui+1 from an initial condition of Ui.

 

[bigduke88]

Hello everyone
It took a long time, but I had to do some other things, about my project. The project describes the process of diffusion of an acid gas, with a freely falling drop of an alkalyne solution. One of the basics in the transport phenomena.
That's why I was talking about dimensionless times etc
Here is source, that I've written about the equation describing the velocity of the falling drop in function of time.
At the end, I've implemented a "Z", this is the distance passed by the drop.
The code has some errors and now I'm trying to fix them. When I do it, will be posted here.
Best regards!!!
P.S.
The source code

IMPLICIT DOUBLE PRECISION (A-H,O-Z)
open(unit=3,file='du/dt.txt')
NIT=200
DT1=1D-8
DT=DT1
DTau=4.D0*DT
Dp=0.00271D0
Ro=998.2D0
Mu=0.001005D0
Diff=2.97271D-9
Tau=(4.D0*Diff*T)/(Dp**2.D0)
T=(Tau*(Dp**2.D0))/(Diff*4.D0)
imax=100
DO 30 i=1,imax
30 U(i)=0.
T=0.
DO i=1,imax
T=T+DT
Tau=Tau+DTau
Re=(Dp*U(i)*Ro)/(Mu)
Cd=(24/Re)*(1+0.125*(Re**0.72))
U(i+1)=U(i)*(1-0.7366*U(i)*Cd*DT)+0.998*DT
DO i=1,imax
U(i+1)=U(i)*(1-0.7366*U(i)*Cd*DT)+0.998*DT
Z=Z+DT*(U(i+1)+U(i))

100 continue
write(18,*) Z,T,Tau
STOP
END

 

[salgerman]

Yeah, I can see a couple of errors...

First, you have not declared the array U(). When you do, you can simply declare it of fix length like U(100), if you are going to stick to imax=100.

Second, you have "100 continue", but the "100" label is missing from your "DO" statements.

I would recommend you stop using numeric labels, altogether, and start using the "DO..END DO" construct. As it is, you have 2 do-loops ended using the same exact continue...not cool, anymore.

Also, in Fortran 90, you can deal with arrays a-la-matlab; in other words, you do no need to initialize U() one item at a time as in

      DO 30 i=1,imax
30      U(i)=0.

Instead, you can simply write

      U = 0.

Lastly, I don't see why you have to re-calculate Cd every time inside the loop if it is a constant.

my 2 cents.

 

[salgerman]

oops...sorry, I guess Cd is changing.

Also, I am confused about your nested DO-loops and you re-calculating the entire U() array in the inner one...without know anything, something does not look right...

 

[gummibaer]

Well,
tracking down bugs is sure a tedious business - but the best way to learn programming.

But one (or two) hint(s) as to good programming practice:
Never use implicit type declaration - unless you want to train your eyesight and attention by tracking down typos. It is more work at first - but I'd recommend a data dictionary anyway. Or are you really sure you still know what you thought today when you have to edit your code in six months ?

Norbert.


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[bigduke88]

Really thank you all for the responses, they are really usefull!!!
I've done some corrections on the source, now it works, but when I open the folder, where should be stored the results from the calculations, there are only F.files and application...
I thought, that when I enter the statement"open(unit=3,file='du/dt.DAT') i could find a .dat file with the results.
The work will continue
Best regards to all!!!

 

[bigduke88]

Excuse me for the repost, but I've forgotten to post the new source code. If it is interessting, I could explain, what is the exact purpose of my Fortran "battle".
The source code:
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION U(100)
open(unit=3,file='du/dt.DAT')
NIT=2000
jprint=500
DT1=1D-8
DT=DT1
DTau=4.D0*DT
Dp=0.00271D0
Ro=998.2D0
Mu=0.001005D0
Diff=2.97271D-9
Tau=(4.D0*Diff*T)/(Dp**2.D0)
T=(Tau*(Dp**2.D0))/(Diff*4.D0)
imax=100
U(i)=0.
T=0.
DO 100 i=1,imax
T=T+DT
Tau=Tau+DTau
Re=(Dp*U(i)*Ro)/(Mu)
Cd=(24/Re)*(1+0.125*(Re**0.72))
U(i+1)=U(i)*(1-0.7366*U(i)*Cd*DT)+0.998*DT
U(i+1)=U(i)*(1-0.7366*U(i)*Cd*DT)+0.998*DT
Z=Z+DT*(U(i+1)+U(i))

100 continue
write(18,*) Z,T,Tau
STOP
END

 

[gummibaer]

The filename you selected seems somewhat inappropriate ('/' is not legal in filenames). I wonder why you do not receive a runtime error on this.

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[salgerman]

Nothing wrong with "/" in file name...just think of it as being in a sub-directory.

bigduke: is the *.dat file in a subfolder named "du" ?

 

[gummibaer]

ohh, that's new to me. I allways used '\' for subdirectories instead of '/'

Learning never ends....

Norbert

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[bigduke88]

Hello everyone
salgerman - yes I've created a subfolder named "dU_dT", because windows does not accept the "/". At this folder I can just find some .dat files, an application file named source and an F file. I was wondering, why there is no text file with the reuslts of the calculations...
This is the newer version of the source code:
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (NR=100)
DIMENSION U(NR)
open(UNIT=13,STATUS='new',file='dU\dT.TXT')
NIT=2000
jprint=500
KPRINT=1
DT1=1D-8
DT=DT1
DTau=4.D0*DT
Dp=0.00271D0
Ro=998.2D0
Mu=0.001005D0
Diff=2.97271D-9
Tau=(4.D0*Diff*T)/(Dp**2.D0)
T=(Tau*(Dp**2.D0))/(Diff*4.D0)
imax=NR
U(i)=0.
T=0.
DO 100 i=1,imax
T=T+DT
Tau=Tau+DTau
Re=(Dp*U(i)*Ro)/(Mu)
Cd=(24/Re)*(1+0.125*(Re**0.72))
U(i+1)=U(i)*(1-0.7366*U(i)*Cd*DT)+0.998*DT
C DISTANCE
Z=Z+DT*(U(I)+U(I+1))

100 continue
write(UNIT=13,FMT=*) Z,T,U
STOP
END


And there are no problems about running it, but there are no results and I don't know why.
The results should be separated in three colums: one for the distance "Z", one for the time "T" and one for the velocity "U".
Now I'm going to work, because there was plenty of other work today.
Best regrads!!!

 

[gummibaer]

You still have the '\' in the open-statement. If '\' creates a subfolder, then your program creates the subfolder. In your current directory you should have a directory called dU and in there you may find a file dt.txt. At least that is what I would expect from your code (with what I have learned ).

Rewrite your open to something like

open (unit = 13, file = 'dU_by_dT.txt')

and after your prog executed you will find a textfile of that name in your current directory. Note, I skipped the 'new' clause, because you will receive a runtime error for every run of your program when you did not delete your outputfile from the run before first. Could be annoying while debugging.

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[salgerman]

bigduke:

Hhhmmm, there might be some kind of confusion here...

I am not sure if YOU mean for "\" to be part of the filename alone...which it wouldn't...I think Windows is interpreting "\" as part of a relative full path name, which includes a subdirectory "dU" and a file named "dT.TXT".

So, I will ask again...do you have a directory named "dU" ?!! Not "dU_dT", but "dU"?

Is the program actually running? Or is it crashing because it cannot find the file "dU\dT.TXT"?

Are you running the program from the command line to see feedback or are you simply double-clicking on the executable?

Also, follow Norbert's suggestion about the name change...that should most probably clear up your problem.

 

[bigduke88]

Thank you, I've done this : open(UNIT=13,STATUS='new',file='dU_by_dT.TXT')
and at the end of the program
write(UNIT=13,FMT=*) 'Z,T,U'
write(unit=13,FMT=*) Z,T,U
And it creates a text document named dU_by_dT, but there are only Z,T,U as head of column, but there are no data. When I start the program, it initiates a window, where normally I colud see the results of the iterations, but this "window" just appears for a moment and then, nothing happens.
Maybe the mistake is in the print statements?

 

[salgerman]

Please post the entire code, once again...the one above has the write statement outside the LOOP!